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ABSTRACT 

Or 



The G333 giant molecular cloud contains a few star clusters and Hn regions, 
plus a number of condensations currently forming stars. We have mapped 13 of these 
sources with the appearance of young stellar objects (YSOs) with the Spitzer Infrared 
Spectrograph in the Short-Low, Short-High, and Long- High modules (5-36 pm). We 
use these spectra plus available photometry and images to characterize the YSOs. 
The spectral energy distributions (SEDs) of all sources peak between 35 and 110 pm, 
thereby showing their young age. The objects are divided into two groups: YSOs 
CN| ■ associated with extended emission in Infrared Array Camera (IRAC) band 2 at 4.5 pm 

('outflow sources') and YSOs that have extended emission in all IRAC bands peaking 
at the longest wavelengths ('red sources'). The two groups of objects have distinctly 
different spectra: All the YSOs associated with outflows show evidence of massive 
envelopes surrounding the protostar because the spectra show deep silicate absorption 
features and absorption by ices at 6.0, 6.8, and 15.2 urn. We identify these YSOs with 
massive envelopes cool enough to contain ice-coated grains as the 'bloated' protostars 
in the models of Hosokawa et al. All spectral maps show ionized forbidden lines and 
polycyclic aromatic hydrocarbon emission features. For four of the red sources, these 
lines are concentrated to the centres of the maps, from which we infer that these YSOs 
are the source of ionizing photons. Both types of objects show evidence of shocks, with 
most of the outflow sources showing a line of neutral sulphur in the outflows and two 
of the red sources showing the more highly excited [Ne in] and [S iv] lines in outflow 
^ • regions at some distance from the YSOs. The 4.5- pm emission seen in the IRAC band 

2 images of the outflow sources is not due to H2 lines, which are too faint in the 5 - 
10 pm wavelength region to be as strong as is needed to account for the IRAC band 
2 emission. 

Key words: stars: formation - stars: massive - stars: protostars - ISM: jets and 
outflows - infrared: ISM - infrared: stars. 



1 INTRODUCTION 

Massive stars form from gas and dust condensations in 
molecular clouds. Although there have been significant ad- 
vances in our understanding of massive star formation, out- 
standing questions remain. No one at this time is able to 
compute the evolution of a star from molecular cloud clump 
all the way to the main sequence because of the huge range 
of scales involved. Numerical models of molecular clouds col- 
lapsing under the force of gravity assume that when enough 
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mass gets into the smallest grid cube or cubes, a star is 
formed. These models, with fragmentation according to the 
local Jeans mass, produce both small and massive stars in 
agreement with observations (e.g., Bonnell & Bate 2006; Pe- 
ters et al. 2010; Wang et al. 2010). 

A serious problem in understanding the formation of 
massive stars had long been that the radiation pressure at 
the surface of any star greater than about 8 Mg should 
prevent the further addition of any mass, and yet stars an 
order of magnitude more massive are known. However, it 
has now been shown through three-dimensional models of 
the collapse of a rotating core that the actual accretion can 
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occur through an optically-thick accretion disc and the ra- 
diation escapes along the axes (e.g., Krumholz et al. 2009), 
thus solving the problem of the excessive radiation pressure 
keeping the infalling material off the star. Computations of 
low-mass stars (Vorobyov & Basu 2005, 2006, 2010) show 
that the disc fragments due to gravitational instabilities and 
the fragments lose angular momentum through viscosity and 
traverse through the disc until very close to the star, at 
which time it is assumed they are accreted. Although such 
episodes of increased mass accretion are observed for low- 
mass stars (e.g., the FU Orionis stars, Hartmann & Kenyon 
2006), such episodic accretion has not been observed for 
massive stars, probably because most mass is accreted while 
the star is still completely obscured optically and it is diffi- 
cult to distinguish the star itself at mid-infrared (MIR) and 
far-infrared (FIR) wavelengths from the warm dust of the 
natal cloud. Other computations indicate that such fragmen- 
tation and resulting accretion are also likely for high-mass 
stars (Kratter, Matzner, & Krumholz 2008). On the other 
hand, the mechanism of accretion onto the surface of the star 
from the disc is very complex and is still not understood (see 
the review of McKee & Ostriker 2007) . 

The internal structure and evolution of an accreting 
massive protostar has been modelled by Hosokawa, Yorke, 
& Omukai (2010), who find that when a model that is ac- 
creting at a rate of 10 -3 Mq yr _1 approaches ~ 10 M Q , it 
'bloats up' to a radius ~ 100 Rq . They predict that proto- 
stars at that stage are unable to ionize any H n region, in 
spite of their high luminosities, and that this stage should be 
observable as a high-luminosity young stellar object (YSO) 
with no Hn region. 

A molecular cloud collapsing into stars is by its nature 
turbulent. It has been suggested that supersonic turbulence 
is crucial in both preventing precipitous collapse on large 
scales and promoting the density enhancements needed at 
smaller scales to actually form stars (e.g., Zinnecker & Yorke 
2007; McKee & Ostriker 2007). Turbulence is induced into 
molecular clouds both on the large scale by rotation through 
a Galactic spiral arm and on the small scale by outflows from 
newly formed stars (e.g., Wang et al. 2010; but see Padoan 
et al. 2009). Given enough time, supernovae are major con- 
tributors. 

The giant molecular cloud (GMC) complex centred at 
/ = 333.2, b = -0.4 (the G333 cloud) is ideal for testing 
theories of star formation. We have been studying the ef- 
fect of turbulence on giant molecular cloud structure and on 
star formation through molecular line mapping of the G333 
cloud (Bains et al. 2006; Wong et al. 2008; Lo et al. 2009). 
This GMC contains the very luminous H II region and clus- 
ter G333.6— 0.2, plus a number of other, smaller star clus- 
ters. The major star clusters, visible optically or in the near- 
infrared (NIR), all lie along the major axis of the cloud. The 
FIR (Karnik et al. 2001) and mm (Mookerjea et al. 2004) 
maps and 6.7 GHz Class II methanol maser surveys (Walsh 
et al. 1998) also locate regions currently forming stars in the 
G333 GMC. 

This paper discusses twelve additional massive star 
forming regions as observed with the Spitzer Space Tele- 
scope (Werner et al. 2004). These are seen by the presence 
of very red, extended sources in images taken for the Spitzer 
Legacy programmes GLIMPSE (Churchwell et al. 2009) and 
MIPSGAL (Carey et al. 2009) with the Spitzer Infrared Ar- 



ray Camera (IRAC, Fazio et al. 2004) and the Multiband 
Imaging Photometer for Spitzer (MIPS, Rieke et al. 2004), 
respectively. Although the usual expectation is that star for- 
mation occurs in the centres of GMCs where the density is 
highest, most of the star-forming regions newly detected in 
the MIR Spitzer images are isolated and lie in the outer 
regions of the cloud. 

Our observed sources were chosen from two groups: (1) 
outflow sources as identified by either the presence of dif- 
fuse emission in the IRAC 4.5 u.m band (these have been 
described as 'extended green objects', hereinafter EGOs, by 
Cyganowski et al. 2008, or 'green fuzzies' by Chambers et 
al. 2009) or the presence of 7 mm SiO emission (Lo et al. 
2007), both thought to arise in shocks; and (2) sources with 
diffuse red emission in the IRAC 8.0-u.m images. All objects 
are compact and bright in the MIPS 24-(im images. All but 
one of the outflow sources also produce 6. 7- GHz Class II 
methanol masers. Fig. 1 shows these objects plotted on the 
IRAC 8.0-u.m image of the G333 cloud (all images in this pa- 
per are log scale). Their coordinates and estimated distances 
(all kinematic assuming the solar distance to the Galactic 
Center, Ro = 8.3 kpc, from Brunthaler et al. 2011 and the 
Galactic rotation curve of Levine, Heiles, & Blitz 2008) are 
given in Table 1 (note that a few sources lie outside the G333 
cloud and appear to be associated with a more distant spiral 
arm) . In this paper we use the word 'sources' to refer to the 
whole extended emission region that we mapped, 'YSOs' as 
the approximately point source including the envelope and 
disc that is the exciting source of the outflow and the MIR 
and FIR emission, and 'protostar' that is the core object 
that is accreting mass and will eventually become a star. 

Our goal for this Spitzer survey is to characterize these 
YSOs according to masses, outflow parameters, and ages. 
We also hope to illuminate observationally some of the the- 
oretical studies mentioned in the previous paragraphs. In 
Section 2 we describe the observations, and in Section 3 we 
describe the results. The spectra of the sources also fall into 
two distinct groups: objects with outflows and ice-absorption 
features at the central YSO, and objects with localized for- 
bidden lines indicating the presence of sources of ionizing 
photons. In Section 3 we also discuss individual objects and 
the frequent multiple sources in a single map, and in Section 
4 we discuss the possible evolutionary stages as indicated by 
the presence or absence of various spectral features. Finally, 
in Section 5 we present our summary and conclusions. 



2 OBSERVATIONS 

We mapped each of our sources with the Infrared Spectro- 
graph (IRS, Houck et al. 2004) with the Short-Low (SL, 5.2 
- 14.5 u.m), Short-High (SH, 10 - 19.5 |am), and Long-High 
(LH, 19.5 - 38 urn) modules on 2009 April 29 - May 7 (a 
week before end of cryogen). The observations were split 
into two groups according to date (2009 April 29-30 and 
May 3-7 May) because the IRS campaign was split by an 
IRAC mini-campaign. Fig. 2 gives an example of the indi- 
vidual IRAC and MIPS images of a source with overlays of 
the IRS slit positions; the exposure identifier number (FITS 
keyword: EXPID) is given for the LH module. Similar fig- 
ures are given for the other sources in the online Supporting 
Information. A background position was observed for each 
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Figure 1. IRAC band 4 (8.0 H-in) image of the G333 cloud in logarithmic scaling. The sources in this survey are marked by square 
boxes. The open circle marks our background position. The outflow sources from the catalog of Cyganowski ct al. (2008) that were not 
observed in our survey are marked by diamonds. Additional red sources that were not observed in our survey arc marked by triangles. 



Table 1. Observational parameters. The table gives for each source the RA and Dec, the date of the observation (end 
of cryogen was 2009 May 15), the AOR identification number, the source type (red or outflow, see text), name of any 
co-located source in the IRAS Point Source Catalog, the kinematic distance estimated from an observed Vlsr., an d 
references: (a) All sources in the G333 cloud are assumed to have the cloud Vlsr = —50 km s" 1 measured by Bains et 
al. (2006). (b) The three sources north of the G333 cloud are assumed to have the Vlsr = ~~ 91 km s~~ 1 measured for 
G333. 168-0.081 by Caswell & Haynes (1987). The masers observed in all these sources have velocities consistent with 
these assumptions (Caswell 2009; Urquhart ct al. 2009). 



Source RA Dec Obs. Date AORKEY Source Type IRAS name Distance Ref. 

(J2000) (J2000) (kpc) 



G332.725- 


-0.620 


16:20:02.8 


-51:00:32 


2009-04-30 


25913088 


Outflow 




3.6 


a 


G332.800- 


-0.595 


16:20:16.4 


-50:56:19 


2009-05-03 


25915392 


Red 




3.6 


a 


G332.813- 


-0.700 


16:20:48.1 


-51:00:15 


2009-04-30 


25912832 


Outflow 


16170-5053 


3.6 


a 


G332.942- 


-0.686 


16:21:18.9 


-50:54:10 


2009-04-30 


25912576 


Outflow 


16175-5046 


3.6 


a 


G332.963- 


-0.679 


16:21:22.9 


-50:52:59 


2009-04-30 


25912064 


Outflow 




3.6 


a 


G332.967- 


-0.683 


16:21:25.0 


-50:53:01 


2009-04-30 


25916160 


Red 




3.6 


a 


G333.131- 


-0.560 


16:21:36.1 


-50:40:49 


2009-05-03 


25916416 


Outflow 




3.6 


a 


G333.163- 


-0.100 


16:19:42.6 


-50:19:53 


2009-04-30 


25915648 


Red 




5.8 


b 


G333.184- 


-0.091 


16:19:45.6 


-50:18:35 


2009-05-03 


25913344 


Outflow 




5.8 


b 


G333.212- 


-0.105 


16:19:56.9 


-50:18:01 


2009-04-30 


25915904 


Red 




5.8 


b 


G333.340- 


-0.127 


16:20:36.9 


-50:13:35 


2009-05-03 


25914624 


Red 


16168-5006 


3.6 


a 


G333.375- 


-0.202 


16:21:06.0 


-50:15:15 


2009-04-30 


25915136 


Red 




3.6 


a 


G333.466- 


-0.164 


16:21:20.2 


-50:09:49 


2009-05-07 


25911552 


Outflow 


16175-5002 


3.6 


a 
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Figure 2. IRAC and MIPS images of G333.131-0.560 and 
G333.128-0.560. Panel (a): IRAC band 1 (3.6 |im) image. Panel 
(b): IRAC band 2 (4.5 (im) image. The two 6.7 GHz methanol 
mascrs (Caswell 2009) arc plotted as black crosses. Panel (c): 
IRAC band 3 (5.8 urn) image with overlays of the SL slit posi- 
tions. Panel (d): IRAC band 4 (8.0 [im) image with overlays of 
the SH slit positions. Panel (e): MIPS 24 nm image with overlays 
of the LH slit positions on which is printed the LH exposure iden- 
tifier numbers. Panel (f): MIPS 70 u.m image. G333.128-0.560 is 
the brightest source in the MIPS 70 um image (panel f) and is 
best located in the MIPS 24 um image (panel e). 



source; the background position is at G333. 3894— 0.9202, or 
RA 16 h 24 m 22!0, Dec. -50° 45' 2". 

The data were reduced and calibrated with the S18.7 
pipeline at the Spitzer Science Center (SSC). The basic cal- 
ibrated data (BCD) from the Spitzer pipeline consist of 128 
by 128 pixel arrays showing 10 spectral orders for the high 
resolution modules and 2 spectral orders for the low res- 
olution modules. For each telescope pointing and for each 
module, we took four spectra (resulting in four BCD im- 
ages) such that the total integration time per pointing per 
module was 24 s. These four BCD images of each detec- 
tor array for each telescope pointing were median-combined 
and the median of all the background BCD images in the ap- 
propriate time group was subtracted. The final backgrounds 
have less noise than any of the source BCD images. This pro- 
cess of subtracting fairly low-noise background BCD images 
greatly reduces the number of 'rogue' pixels 1 (pixels with 



temporary extra dark current or non-linear response as a 
result of hits by cosmic rays or solar protons). The resulting 
background-subtracted BCD images were further cleaned of 
rogue pixels using the interactive version of irsclean 2 . The 
noisy order edges (the ends of the slits) were cleaned with 
our own software. 

Two tools employing the Interactive Data Language 
(idl) from the SSC were used to obtain spectra from the 
BCD images: the Spectroscopy Modeling, Analysis and Re- 
duction Tool (smart, Higdon et al. 2004) and the CUbe 
Builder for IRS Spectral Mapping (cubism, Smith et al. 
2007a). 



2.1 Spectra and line fluxes — SMART 

SMART can be used to extract spectra from individual BCD 
images, or sections of BCD images for long slit spectra. We 
used SMART to extract spectra from the SH and LH BCD im- 
ages. Because the IRS was calibrated using point sources but 
all our sources are extended, a correction for the telescope 
diffraction pattern as a function of wavelength is necessary; 
these 'slitloss' correction factors were taken from tables pro- 
vided by the SSC . SMART returns separate spectra for each 
of the ten SH or LH spectral orders. Many of the spectra in 
certain orders, particularly LH, are contaminated by spec- 
tral fringing of order 3-6 per cent. An advantage of the use 
of smart to measure line fluxes is that it includes a tool 
(irsfringe) to fit and subtract the fringes, order by order. 

Line and continuum fluxes were estimated by fitting 
Gaussian profiles plus a linear continuum to each of the 
line profiles in the spectra using both idl's CURVEFIT pro- 
cedure and the procedures in SMART. For the faintest lines, 
the fluxes were estimated by integrating over the line profile. 
Errors for all fluxes were estimated by computing the root- 
mean-squared deviation of the data from the fit; because of 
the systematic effects of low-level, uncorrectable rogue pix- 
els on the spectra, these errors are consistently larger than 
errors estimated any other way, such as the statistics of me- 
dian averaging or the counting statistics plus read noise from 
the SSC's pipeline. 

The lines discussed in this paper are listed in Table 2, 
along with references to the atomic or molecular physics 
parameters (wavelengths, transition probabilities, collisional 
cross sections) needed for further analysis. 



2.2 Spectra and line fluxes - CUBISM 

CUBISM (Smith et al. 2007a) combines all the spectra from 
a single or multiple Astronomical Observation Requests 
(AORs) for a given module if SH or LH or a given order 
if SL (SL order 1 covers 7.5 - 14.2 urn and order 2 covers 
5.2 - 7.5 u.m) or LL (Long-Low, 14 - 38 (im) and produces a 
three-dimensional cube of x, y, and wavelength, where x and 
y refer to coordinates parallel and perpendicular to the mod- 
ule slit. It was designed for spectral maps of galaxies from 
the SIRTF Nearby Galaxies Survey (SINGS, Kennicutt et 
al. 2003) and the calibration assumes that the source is ex- 
tended. The user selects an area of the sky in x and y to 



1 http://irsa.ipac.caltech.edu/data/SPITZER/docs/irs/features/ 



2 http://irsa.ipac.caltech.edu/data/SPITZER/docs/dataanalysistools/tools/irsclca 

3 http://irsa.ipac.caltech.edu/data/SPITZER/docs/dataanalysistools/tools/spice/ 
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Table 2. Line parameters. For each line the table gives the wavelength, the ionization poten- 
tial (IP, the energy required to produce the ground state of the molecule or ion), the energy 
of the upper level of the transition, the extinction law for that wavelength as a fraction of 
the extinction law for 9.6 urn, and relevant references. The references in the last column are 
as follows: (1) Bautista & Pradhan (1996), (2) Dufton & Kingston (1991), (3) Erickson ct 
al. (1989), (4) Feuchtgruber ct al. (1997), (5) Griffin, Mitnik, & Badnell (2001), (6) Kelly 
& Lacy (1995), (7) McLaughlin & Bell (2000), (8) Nussbaumcr & Storey (1988), (9) Pelen 
& Berrington (1995), (10) Quinct (1996), (11) Quinct, Lc Dorneuf, & Zeippen (1996), (12) 
Ramsbottom et al. (2005), (13) Saraph & Storey (1999), (14) Storey & Hummer (1995), (15) 
Tayal (2004), (16) Tayal (2006), (17) Tayal & Gupta (1999), (18) Wolnicwicz, Simbotin, & 
Dalgarno (1998), (19) Zhang (1996). 



Line 


Wavelength 


IP 


Eupper 


Tx/rg.6 


References 




(urn) 


(eV) 


(cm- 1 ) 






H 2 S(0) 


28.219 





354.37 


0.369 


18 


H 2 S(l) 


17.035 





705.69 


0.502 


18 


H 2 S(2) 


12.279 





1168.80 


0.344 


18 


H 2 S(3) 


9.665 





1740.36 


0.997 


18 


H 2 S(5) 


6.909 





3187.64 


0.273 


18 


H 2 S(7) 


5.512 





5002.04 


0.290 


18 


Hi 7-6 


12.372 


13.60 


107440.45 


0.332 


14 


[Oiv] 2 P 3/2 - 2 Pi/ 2 


25.890 


54.94 


386.245 


0.404 


4, 16 


[Neil] 2 P 1/2 - 2 P 3/2 


12.814 


21.56 


780.42 


0.313 


5 


[Nera] 3 Pi- 3 P 2 


15.555 


40.96 


642.88 


0.417 


4, 7 


[Sin] 2 P 3/2 - 2 Pi/ 2 


34.815 


8.15 


287.23 


0.306 


2, 4 


[Siii] 3 Pi- 3 P 


33.481 


23.34 


298.68 


0.311 


4, 17 


[Siii] 3 P 2 - 3 Pi 


18.713 


23.34 


833.06 


0.548 


6, 17 


[Siv] 2 P 3/2 - 2 P 1/2 


10.511 


34.79 


951.43 


0.777 


13 


[Am] 2 P 1/2 - 2 P 3/2 


6.985 


15.76 


1431.58 


0.272 


9 


[Fen] 6 D 7/2 - 6 D 9/2 


25.988 


7.90 


384.79 


0.403 


1, 8, 11, 12 


[Fein] 5 D 3 - 5 D 4 


22.925 


16.19 


436.21 


0.448 


3, 10, 19 



extract a spectrum. Since the output spectrum header con- 
tains the coordinates on the sky of the corners of the x,y 
rectangle, there is an option to extract a spectrum from a 
different module of exactly the same area of the sky. 

We used CUBISM to first extract the spectrum of every 
telescope pointing observed by the SH module (e.g., panel 
'd' of Fig. 2) and then extracted spectra from the LH and 
SL cubes of all the same positions on the sky (skipping those 
map corners that were not observed by all three modules). 
To increase the signal/noise ratio (S/N), we also extracted 
spectra from regions that looked especially interesting, such 
as large areas seen to have H II region lines, appropriate ar- 
eas for background subtraction, and point sources where the 
observed grid did not exactly overlay the YSO. Because the 
spectra covered exactly the same regions of the sky in all 
modules, they could be easily stitched together to make a 
complete 5.2 - 35 spectrum with very little flux adjust- 
ment at the wavelengths where the modules overlap (e.g., 
Smith et al. 2007b). 

The spectra of the YSOs are shown in Fig. 3. These 
spectra will be discussed further in Section 3, but we wish 
to point out the following features: All spectra are domi- 
nated by a warm continuum that increases steeply to longer 
wavelengths. The sharp lines are the Hn region and pho- 
todissociation region (PDR) lines listed in Table 2; these 
are present in almost every spectrum, including the 'back- 
ground' spectrum, and are probably due to the same diffuse 
interstellar medium (ISM) that is seen in the molecular cloud 
image of Fig. 1. The broad emission features are interstellar 
polycyclic aromatic hydrocarbon molecules (PAHs) at 6.2, 



7.7, 8.6, 11.2, 12.0, and 17 um plus additional wavelengths 
as indicated in the figure; they are also seen in almost all 
spectra, including the background. The absorption features 
include silicates at 9.6 and 18 um and ices at 6.0 (H 2 and 
HCOOH), 6.8 (CH 3 OH and NH+) and 15.2 um (C0 2 ) (e.g., 
Boogert et al. 2008). 

The advantage of CUBISM over SMART is that a complete 
spectrum can be produced for a given region of the sky and 
there are no 'jumps' in the spectrum at the SH and LH 
order boundaries (which jumps are due to smart's use of 
point sources for calibration) . The advantages of SMART over 
CUBISM is that it can be run as a script for all SH and LH 
spectra and the LH spectrum that was extracted by SMART 
can be defringed whereas the 19 - 35 u.m LH spectrum that 
was extracted by CUBISM can not. We will give an example 
of the need for defringing in Section 3.6. 



3 RESULTS 

We start with the basic results for all sources (extinction, 
luminosity) and then detail the results from the ice features 
and the forbidden line analysis. We end with a brief discus- 
sion of each YSO. 



3.1 Extinction - PAHFIT 

The MIR extinction was determined by fitting the observed 
spectra with the program pahfit, developed by Smith et 
al. (2007b) and available online. We use the fitted 9.6 \im 
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Figure 3. Panel (a): spectra of 'red' sources with identification 
of features of interest. Panel (b): spectra of 'outflow' sources with 
identification of features of interest. 



optical depth and MIR extinction law to correct the mea- 
sured line intensities, both hydrogen and forbidden lines, 
for extinction and to determine the variations in extinction 
between the sources. 

The extinction in the MIR (> 8 \ixn) is due to absorp- 
tion by silicates. This is most notable in the depth of the 
9.6 um silicate feature in some of the spectra in Fig. 3, al- 
though in many of the spectra this feature is obscured by the 
gap in the PAH spectrum from 9-11 |j.m. Consequently, the 
whole spectrum needs to be fitted at the same time: MIR 
continuum, PAHs, and the 9.6 and 18 |J.m silicate absorption 
features. 

Smith et al. (2007b) developed a program, called pah- 
fit, to fit all these components by non-linear least-squares to 
the low-resolution SINGS galaxy spectra. This idl program, 
with modifications, was used to estimate the MIR extinction 
affecting all our spectra. Although pahfit could produce 
good fits for all our spectra that are most like H n region 
spectra (the red sources and the outskirts of all the maps, 
e.g., Fig. 3a), as might be expected since the SINGS spiral 
galaxy spectra are dominated by H n regions, the spectra of 
the YSOs (e.g., Fig. 3b) that are dominated by absorption 
features often gave poor fits, pahfit uses a Drude profile 
to fit each of the 24 PAH features in its database; however, 
Drude profiles have very broad wings, and the least-squares 
fitting routine sometimes fitted the 5-10 [im continuum 
with exceptionally strong PAH features and no 5 - 10 U-m 
continuum. This was particularly likely for the sources with 
the 6.0 and 6.8 urn ice features. (A non-physical result is 
something all users of least-squares fitting routines must al- 
ways watch out for.) 

Because most of the PAHs in our spectra arise in the 
diffuse ISM, we expect, in fact, that the PAH contribution 
should be similar for all the spectra (there can be substantial 
variations in the PAH spectrum between types of sources; 
see e.g., the review of Tielens 2008). Consequently, in or- 
der to improve the fit and to fit all sources consistently, we 
modified the pahfit program to use a template for the PAH 
features instead of the Drude profiles for each individual fea- 
ture. (We kept the structure of the pahfit program rather 
than write a new program because of all its other useful fea- 
tures: silicate absorption, continuum fits, line fits, graphics, 
etc.) Although it would be preferable to use a PAH tem- 
plate from the G333 region, none is available that does not 
include at least some extinction, since the G333 cloud is at 
a distance ~ 3.6 kpc (Table 1). An unextincted PAH spec- 
trum of high S/N that has minimal forbidden lines taken 
with Spitzer's IRS SL module is needed for the template - 
we found such a spectrum in Spitzer program 0045 (PI: T. 
Roellig) just south of the Orion Nebula bar at 5 h 35 m 2r!3 
— 5°25'6". We extracted this spectrum with CUBISM over a 
72-arcsec 2 aperture from SL and SH AORs 4117760 and 
4118016, stitched on a longer wavelength contribution from 
the Midcourse Space Experiment (MSX) spectrum of the 
whole Orion Nebula (Simpson et al. 1998), and estimated 
the continuum with pahfit. Fig. 4 shows the spectrum and 
fit; the PAH template consists of the data (square boxes) 
minus the fitted continuum (solid black line) over intervals 
between 5.2 and 17.6 \im (see Smith et al. 2007b for a de- 
scription of the other fit components). This fit found zero 
silicate absorption. 

For the MIR extinction law, we used the Galactic Centre 
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Figure 4. Spectrum of the Orion Bar from Spitzcr program 0045 
with results from PAHFIT overplotted. The data, marked by square 
boxes with error bars, consist of IRS SL spectra from 5.24 to 
12.44 urn, SH spectra smoothed to the SL resolution from 12.44 
to 19.05 |xm, and spectra from the full Orion Nebula measured 
by MSX (Simpson et al. 1998) from 19.05 to 25.65 urn. Lines of 
[Neil] 12.8 urn, [Cin] 14.4 urn, [Neni] 15.6 urn, and [S III] 18.7 
were excised from the spectrum before fitting. The results from 
PAHFIT are as follows: the fit is plotted in green, narrow lines are 
violet, PAH features are blue, blackbody continua are red, and the 
total continuum is plotted in black. The PAH template consists 
of the Orion data minus the total continuum from 5.24 - 8.85 urn, 
10.77 - 14.41 urn, and 16.29 - 17.59 urn; from 8.85 - 10.77 urn 
and 14.41 - 16.29 urn the template consists of the wings of the 
PAHFIT PAH features, which present a much less noisy continuum 
than the Orion data. 
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Figure 5. Extinction law. The solid line is the extinction law 
as was used for the G333 cloud; it was taken originally from the 
Galactic Centre law of Chiar & Tielens (2006) but modified in 
the 13.6 - 16.8 um region (see text) and at wavelengths short of 
8 |xm where we use the extinction law of Indebetouw et al. (2005). 
The original law of Chiar & Tielens is plotted as the dashed line. 



law of Chiar & Tielens (2006), slightly modified to include 
a little more extinction at 15 urn. This is shown in Fig. 5, 
normalized to 1.0 at 9.6 urn. The reason for the modifica- 
tion is that otherwise, our very heavily extincted positions 
would have a little bump in the fitted spectrum at 15 urn 
that is not observed. The modification flattened the bump; 
however, a larger increase in the extinction might actually 
be the case but the exact shape of the extinction curve at 
this wavelength would be difficult to determine because the 
objects with the largest extinction also exhibit the 15.2 urn 
C0 2 ice feature. The Chiar & Tielens (2006) Galactic Centre 
extinction has a larger 18/9.6 um ratio than the extinction 
laws of Draine (2003a,b); we prefer it because it gives a bet- 
ter fit to the H2 S(0), S(l), and S(2) line ratios, which lines 
arise in the diffuse ISM along the line of sight to the Galactic 
Centre (Simpson et al. 2007). 

The final modification was to add a component to the 
fitting procedure for the ice features. To keep to the struc- 
ture of the pahfit program, the ice features were added as 
Gaussian absorption lines at 5.93, 6.80, and 15.30 urn with 
full width at half maximum (FWHM) of 0.429, 0.843, and 
0.602 urn. Gaussian line profiles are certainly not appropri- 
ate, but all that is needed for this modification is to make 
some accounting for the strength of the absorption relative 
to the continuum and the PAH emission, since the purpose 
of our using pahfit is to estimate the extinction and the 
strength of the PAHs. The feature central wavelengths and 
FWHM listed above were estimated using the line fitting 
procedure in SMART. 

The pahfit deconvolutions are given in Figs. 6 and 7 
for the spectra at the YSO positions of all the sources. The 
optical depths at 9.6 urn, Tg.e, are at the YSO positions. 
Note that compared to the spectra in Fig. 3, we used the SL 
spectra for the 5.2-14 urn wavelength region and degraded 
the resolution of the rest of the SH and LH spectra when 
we stitched all spectra together in order to better match 
the line widths in pahfit, which assumes all spectra are 
from the low resolution modules. We remark here that the 
PAH template is never a perfect fit to the observed spectra, 
and that each source deviates in a different way, regarding 
both feature strengths and shapes. See the review of Tielens 
(2008) for more discussion of the possible components to the 
PAH spectrum. 

The computed values of T9.6 plotted as functions of po- 
sition are given in Figs. 8 to 18 for all maps except G333.131 
(Fig. 2), which has too little S/N for production of stitched 
spectra from all three modules. The values of T9.6 are also 
given in Table 3 for each source, both the value of rg.6 at the 
central YSO (this is usually the deepest absorption) and a 
value of rg.6 from the outskirts of the spectral map for an es- 
timate of the background extinction. The sources with the 
largest distances from the Sun (Table 1) have the largest 
background extinctions as might be expected if the back- 
ground extinction is due to the diffuse ISM. We also give the 
difference in extinction of source minus background. Here 
we see that for the most part the outflow YSOs have larger 
extinction than the red sources. We note that we assumed 
in the pahfit computation that the unextincted continuum 
consists of a smooth black body with no silicate emission. 
This should be appropriate for a very optically thick YSO 
core, but it is probably not correct for the diffuse ISM of 
the background positions. Including silicate emission as well 
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Figure 6. Fit results for the SL and smoothed SH data for the po- 
sition that includes the YSO for the red sources plus the G333.466 
position that has the maximum forbidden line emission. The data 
are plotted in black, the PAH template, as absorbed by the extinc- 
tion curve multiplied by the rg.6 as given in each panel, is plotted 
in blue (dashed), the total fit is plotted in green (see text). The 
data minus the extincted PAH template are plotted in red. Note 
that PAH template is never a really good fit, but the deviations 
vary with the different sources, indicating that there is no one 
single PAH spectrum, even for very similar sources. Panel (g): 
G333.466 is an outflow source that has an Hn region component 
to the west of the YSO and outflow. The plotted position con- 
tains the peak of the Hn region component, and shows a typical 
Hll region spectrum as seen in panels (a) to (f), except that it 
also has a substantial CO2 15- urn ice absorption feature from the 
accreting dust and gas. 



Figure 7. Fit results for the SL and smoothed SH data for the po- 
sition that includes the YSO for the outflow sources. The colours 
are described in Fig. 6. The ice features were represented by Gaus- 
sians in the fit process for lack of a good template and thus their 
fits are not accurate. Panel (d): note that even though the fit- 
ting program found no PAH components for G332.963, there are 
small bumps at 6.2, 7.7, and ff.3 urn in the spectrum, showing 
that some PAH emission (probably foreground) is present, even 
here. 

as silicate absorption could produce significantly larger es- 
timates for the background extinction (Simpson & Rubin 
1990). 

3.2 Luminosities and associated parameters 

The source luminosities can be estimated from the total 
spectral energy distribution (SED) and the distance. How- 
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Figure 8. Extinction map and SH slit positions for 
G332.725-0.620 plotted on the IRAC band 2 (4.5 (im) image. 
Positions without a value for the optical depth do not have SL 
spectra. 



Figure 10. Extinction map and SH slit positions for 
G332. 813-0.700 plotted on the IRAC band 2 (4.5 urn) image. 
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Figure 9. Extinction map and SH slit positions for 
G332.800-0.595 plotted on the IRAC band 2 (4.5 urn) image. 



ever, the luminosity cannot be estimated by simply inte- 
grating over the observed fluxes as a function of wavelength 
because YSOs are not spherically symmetric. Robitaille et 
al. (2006) computed a 'grid' of 20,000 YSO models, each cal- 
culated with 10 inclination angles ranging from cos i = 0.05 
to cos i = 0.95 (i = 87° to i = 18°). The models consisted 
of a protostar, a disc, and an envelope with cavities along 
the axes due to outflows. They found that the optical depth 
to the central protostar can be quite variable depending on 
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Figure 11. Extinction map and SH slit positions for 
G332. 942-0.686 plotted on the IRAC band 2 (4.5 urn) image. 



the axis inclination, i. For pole-on views (i = 0), the central 
protostar would be visible and there can be substantial ob- 
served NIR and even visible flux, whereas for i — 90, the star 
can be completely obscured by an edge-on, optically thick 
disc or envelope (see also Whitney et al. 2004). 

Consequently, to estimate the luminosities we used the 
SED fitter of Robitaille et al. (2007) 4 . This online appli- 
cation takes input fluxes as a function of wavelength, esti- 
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Figure 12. Extinction map and SH slit positions for 
G332.963-0.679 and G332.967~0.683 plotted on the IRAC band 
2 (4.5 |xm) image. 
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Figure 13. Extinction map and SH slit positions for 
G333.163-0.100 plotted on the IRAC band 2 (4.5 urn) image. 

mated distances (Table 1), and visible extinction and fits 
models from the Robitaille et al. (2006) grid. For the esti- 
mated extinction, we used either the estimated background 
extinction from Table 3 or we used a very large value so that 
the extinction is a completely free parameter. As would be 
expected, the fits with the constrained extinction have much 
larger x 2 ; they also exhibit a larger range of fitted parame- 
ters, particularly the inclination angle and the model age. 

The ranges of the fitted parameters are listed in Ta- 
ble 3 for luminosity, L/Lq, protostar mass, M*/Mq, enve- 
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Figure 14. Extinction map and SH slit positions for 
G333. 184-0.091 plotted on the IRAC band 2 (4.5 urn) image. 




60 58 56 54 

RA (J2000) 16 h 19 m 



Figure 15. Extinction map and SH slit positions for 
G333.212-0.105 plotted on the IRAC band 2 (4.5 urn) image. 

lope mass, Menv/M©, protostar effective temperature, T c g, 
inclination, and model age for the best-fitting models. The 
parameters in Table 3 are combined from fits made with con- 
strained and unconstrained extinction. From the fits, which 
are shown in Fig. 19 for the fits with the unconstrained ex- 
tinction, we also estimate and give in Table 3 the wavelength 
of the maximum of the fitted SED as plotted in XF\ space. 

Other parameters from the fitted models are of less use 
in characterizing the YSOs: The models with the larger ages 
usually have discs but the disc may or may not contribute 
substantial flux to the SED at the shortest wavelengths; the 
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Figure 19. Spectral energy distributions for all sources with fits from the online SED fitter of Robitaille et al. (2007). The plotted points 
are our data from 5 to 35 urn, fluxes measured from the IRAC and MIPS, and fluxes taken from the literature (see text). The multiple 
solid lines are the SEDs of the top ten fitting models with unconstrained extinction from Robitaille et al. (2006). (A figure showing the 
models with constrained extinction is very similar.) The dashed line is the SED of the central source of the best-fitting model including 
interstellar extinction but not the extinction from the envelope and disc. There are two SEDs for G333.466: the whole source including 
the fluxes of IRAS 16175-5002 and the SED for the YSO with the ice absorption features. 



models with the shortest ages and lowest T e g (the reddest 
SEDs) are the least likely to have discs. The total visible 
extinction, Av, for the SED is the sum of the circumstellar 
and interstellar values - the largest unconstrained interstel- 
lar extinction is Av ~ 102 (in G333.466 YSO), for which 
the model circumstellar Av ~ 463. We constrained the in- 
terstellar Ay to be < 50 for the YSOs with distance 5.8 kpc 



(this is effectively unconstrained) and Ay < 10 to 30 for 
the YSOs with distance 3.6 kpc (Table 1). In general, for 
a given SED the model circumstellar Av varies much more 
widely than the interstellar Av, even by over two orders of 
magnitude for some of the red sources. 

The model inclination is also not well-determined. For 
the fits with unconstrained extinction, all the outflow YSOs 
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Table 3. Results from fitting models to the spectra and to the SEDs. The optical depth, rg.e, is given for the spectrum of the 
source maximum (rg.6 Src), minimum values near the corners of the map which should represent the diffuse ISM or background 
(rg.6 Bkgr), and the difference of source maximum minus background (719.6 Src— Bkgr). The values of model luminosity L/Lq, 
protostar mass M./Mq, envelope mass M env /MQ, protostar effective temperature T e g, inclination (Inch), and model age are the 
range of the best fitting models from the online SED model fitter of Robitaillc ct al. (2007) (see text). A(SEDmax) is the wavelength 
at which the FIR SED is maximum in Fig. 19. Exponents are written as follows: 2-5c3 equals 2 X 10 3 — 5 X 10 3 . 
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^ Luminosity may be underestimated because the map coverage is incomplete. 
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Figure 16. Extinction map and SH slit positions for 
G333. 340-0.127 plotted on the IRAC band 2 (4.5 um) image. 



Figure 17. Extinction map and SH slit positions for 
G333. 375-0.202 plotted on the IRAC band 2 (4.5 um) image. 



have i = 18°, the minimum of the limited number (10) of 
inclination angles available to the SED fitter. For these fits, 
the best fits (lowest x 2 ) have large interstellar Av\ how- 
ever, since all the outflow YSOs are readily visible at the 



shorter wavelengths (IRAC bands 2 and 3), the YSOs must 
be observed almost pole-on where the envelope extinction is 
smallest. On the other hand, when the interstellar extinction 
is constrained to be smaller (Av < 25), the observed depth 
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Figure 18. Extinction map and SH slit positions for 
G333.466-0.164 plotted on the IRAC band 2 (4.5 urn) image. 



of the 10 [im silicate feature and the small NIR fluxes force 
the fit to have a large optical depth to the protostar via the 
envelope and disc (hence a larger inclination angle) - but 
then the YSO must have a hotter star and a larger model 
age to have the hot disc that produces the NIR flux. Because 
the background extinction could be significantly larger (Sec- 
tion 3.1), we regard the differences between fits using the 
constrained and unconstrained extinction as not significant. 
As emphasized by Robitaille et al. (2007), the fits are not 
unique but serve to demonstrate the range of possible pa- 
rameters. These parameters can almost certainly be better 
constrained by detailed modelling (e.g., Sewilo et al. 2010). 

Because the range of fitted parameters is greatly de- 
creased as more wavelengths are included (Robitaille et al. 
2007), we increased the fitted wavelength range as follows: 
No source is optically visible and only G332.800 has useable 
2MASS (Skrutskie et al. 2006) J, H, and Ks magnitudes. 
Some sources appear to have extended emission at Ks, par- 
ticularly the outflow regions. Since this is probably scattered 
light and the YSO itself either is not or is barely visible, we 
do not use these wavelengths. We measured fluxes from the 
IRAC band 1 (3.6 urn) and band 2 (4.5 \im) images (only 
about half the sources have GLIMPSE catalog band 1 or 
band 2 fluxes) and used our own spectra instead of IRAC 
bands 3 and 4 (5.8 and 8.0 um, respectively). Since our spec- 
tra are greatly contaminated by PAHs (and so would also 
be the photometry in IRAC bands 3 and 4), we first sub- 
tracted the PAH templates fitted by pahfit as described in 
the previous section. This enabled us to measure much more 
reliable fluxes at 5.48, 7.99, 9.59, 12.02, 14.02, 18.02, 20.01, 
25.03, 30.00, and 35.02 um for all YSOs except G332.725, 
for which we were not able to obtain a SL order 1 spectrum 
(7.5 - 14 am) and so substituted 7.40 and 10.18 \ixn fluxes 
for the 7.99 and 9.59 urn fluxes. These wavelengths miss 
the ice absorption features and most of the PAHs except at 
7.99 urn, but the template fit and consequent subtraction 



is not bad at that wavelength. Avoiding these emission and 
absorption features is important because they are not mod- 
elled by Robitaille et al. (2006). Examples of spectra with 
the PAH template subtracted are given in Figs. 6 and 7 and 
the resulting continuum fluxes are given in Table Al. 

For longer wavelengths we measured the flux at 70 um 
on the MIPS images and used Infrared Astronomical Satel- 
lite (IRAS) fluxes where available (Table 1). We used fluxes 
from the AKARI (Murakami et al. 2007) point source cat- 
alog (Yamamura et al. 2010) at 65, 90, 140, and 160 um, 
fluxes from Karnik et al. (2001) at 150 and 210 urn, and 
fluxes from Mookerjea et al. (2004) at 1.2 mm. In the future 
we will add the 150 um to 600 um fluxes estimated from 
the images taken by the Herschel Space Observatory (Pil- 
bratt et al. 2010) as part of the Hi-GAL Open Time Key 
Project (Molinari et al. 2010). These images were taken with 
the Photodetector Array Camera and Spectrometer (PACS) 
(Poglitsch et al. 2010) and the Spectral and Photometric 
Imaging Receiver (SPIRE) (Griffin et al. 2010) cameras; 
they are available online as of early 2011 reduced to 'level 2' 
by the Herschel pipeline but they need much more data re- 
duction for bad pixels and other artefacts before they can be 
used to estimate fluxes. It is beyond the scope of this paper 
to reduce Herschel data, but we can say that all our objects 
can be detected in the Herschel images, including the SPIRE 
600 um images (except for G332.967, which is confused with 
G332.963), and thus the long wavelength parts of the SEDs 
seen in Fig. 19 are reasonable and the YSO envelopes must 
be large as the models predict (Table 3). 

Unfortunately, the longest wavelength fluxes may in- 
clude multiple sources (massive YSOs are, after all, born 
in clusters). The result is that including longer wavelengths 
increases the wavelength of the SED maximum and also in- 
creases the total luminosity. In some cases it may be better 
to regard the luminosities in Table 3 as possibly pertaining 
to a cluster rather than an individual YSO. For example, in 
Table 3 we separate the YSO G333.466 from the total cluster 
IRAS 16175-5002, which has a much larger luminosity. 

3.3 Outflow YSOs with ice absorption features 

All YSOs defined by either their appearance on IRAC 3- 
colour images or by the presence of SiO emission to be out- 
flow sources (Table 1) have ice absorption features at 6.0, 
6.8, and 15.2 |j.m. These spectra are plotted in Fig. 3(b) 
and spectra with the PAH template subtracted in Fig. 7. 
A number of icy compounds have been suggested as pro- 
ducing these features, such as HCOOH and H2O (6.0 um), 
CH3OH and NH+ (6.8 um) and C0 2 (15.2 |im). A review 
of the history of the detection of interstellar ices is beyond 
the scope of this paper - for starters we refer the reader 
to the review by Boogert & Ehrcnfrcund (2004) for a dis- 
cussion of the high resolution spectra taken with the short- 
wavelength spectrometer on the Infrared Space Observatory. 
A great number of spectra of YSOs have been taken with 
Spitzer, see for example, Boogert et al. (2008), Zasowski et 
al. (2009), and Oliveira et al. (2009). With the high sensi- 
tivity of Spitzer, these authors have been able to show that 
there are many absorption features in addition to the main 
features at 6.0, 6.8, and 15.2 (i.m. These new features enable 
the authors to identify additional molecules such as CH4 
at 7.67 um (Oberg et al. 2008) and HCOOH at 7.25 um 
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(Boogert et al. 2008), thus showing details of the chemistry 
of the envelopes and discs of the YSOs. There is a little dip 
at 7.67 \im in most of the outflow sources that may be CHU, 
but our spectra are so overladen with foreground PAHs and 
silicate absorption that it is not possible to estimate its re- 
ality or strength. 

However, we do notice that the fit to the 9.6 \im sili- 
cate feature is never good (Fig. 7) for the outflow sources 
whereas it is quite good for the red sources (Fig. 6). The dif- 
ference can be explained by additional absorption at around 
9.0 urn in the outflow source spectra; an absorption feature 
at this wavelength has been attributed to NH3 ice by Gibb 
et al. (2000). Gibb et al. (2000), Boogert et al. (2008) and 
Bottinelli et al. (2010) also find absorption features due to 
CH3OH ice at 9.7 u.m. Unfortunately, our spectra have so 
much silicate extinction that the signal/noise is very poor at 
this wavelength and thus we cannot detect any such 9.7 |im 
additional absorption. 

We also do not detect the absorption features due to 
gaseous C2H2, HCN, and CO2 that are found in massive 
YSOs in the Galactic Centre by An et al. (2009). Because 
these features are narrow (and hence not confused by PAHs), 
we should be able to detect them if they exist at similar 
strengths to the features measured by An et al. (2009). We 
conclude our objects do not have amounts of warm gas simi- 
lar to the Galactic Centre YSOs studied by An et al. (2009). 
Our objects could be cooler or at an earlier evolutionary 
stage. However, the three YSOs observed by An et al. (2009) 
could still be outflow sources since extended 4.5 u.m emis- 
sion is detected in their vicinities (Chambers, Yusef-Zadeh, 
& Roberts 2011). 

3.3.1 The 15 urn C0 2 ice feature 

It is much easier to reliably detect the CO2 ice feature than 
the features at 6.0 and 6.8 u.m because the 15 (xm region is 
free of the strong PAH emission features seen at the shorter 
wavelengths and at 16.3 - 17.5 [im. As a result, the 6.0 and 
6.8 |J.m features are detected only at or near those YSOs that 
are IRAC band 3 point sources and have substantial 5.5 u.m 
continuum (Fig. 3b). The 15 urn CO2 feature, however, is 
detected more widely, both in the extended region around 
the outflow YSOs and in two red sources (G333.163 and 
G333.375) that have no other ice features. 

Because the CO2 feature can have such good sig- 
nal/noise and is so readily detectable, we have attempted 
to decompose its shape into polar and apolar CO2 ice com- 
ponents as was described by Gerakines et al. (1999). The 
five components of the fit consist of the laboratory mea- 
surements by Ehrenfreund et al. (1997) for ices with ratios 
CO 2 :H 2 O=14:100, CO 2 :CO=26:100, CO 2 :CO=70:100, and 
pure CO2 ice at 10 K. Since there is always smooth absorp- 
tion in our spectra at shorter wavelengths than seen in any 
of Ehrenfreund et al.'s (1997) laboratory spectra, for the 
fifth component we shifted the CO2:CO=70:100 component 
to shorter wavelengths by an arbitrary amount (calculated 
by the least-squares routine). There is no indication of any 
additional 15.4 urn 'shoulder' component as was fitted by 
Pontoppidan et al. (2008) and An et al. (2009), although it 
must be said that our continua are poorly estimated at the 
longer wavelengths of the CO2 feature owing to the PAH 
emission from 16.3 - 17.5 urn. We estimated the compo- 
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Figure 20. Spectra, corrected for extinction, of the six YSOs 
that clearly show the CO2 ice feature with fits to laboratory ice 
spectra from Ehrenfreund et al. (1997). The dashed line is polar 
ice with C02:H20 = 14:100. The other ice features are all apolar 
ice: the dot-dashed line is CC>2:CO = 26:100, the dotted line is 
C02:CO = 70:100, the dot-dot-dashed line is pure CO2 ice, and 
the long dashed line is the CC>2:CO = 70:100 ice feature shifted 
to shorter wavelengths by an arbitrary (calculated) amount (see 
text). The light solid line is the sum of all components and the 
heavy solid line is the fitted continuum, which stops short of the 
16.4 um PAH feature. Only G333.466 (panel f) has any non-zero 
pure CO2 ice component, and there it is small and not statistically 
significant. 



nent strengths by fitting the lab spectra of these five com- 
ponents plus a sloped linear continuum from 14.7 to 16.3 \irci 
by non- linear least-squares using mpfit (Markwardt 2009). 
The polar CO2:H2O=14:100 component at ~ 15.3 u.m is by 
far the dominant, but there is always (if the S/N is ade- 
quate) a shorter wavelength component at ~ 15.1 urn that 
can be attributed to one or all of the apolar components. 
The spectra of the six YSOs with the best S/N are plotted 
in Fig. 20. We list the C0 2 optical depths for these YSOs 
in Table 4, where we have summed the apolar components, 
since there probably is no real physical distinction between 
them. 

The pure CO2 component has substructure with a fea- 
ture at ~ 15.28 urn; that component is identified by the 
least-squares routine only in those spectra that are noisy at 
that wavelength. None of our outflow YSOs reliably show 
this structure in the 15 \im feature, which structure can 
be ascribed to heat-processing of the CO2 ice. Since such 
structure has been seen in certain of the low mass YSOs by 
Boogert et al. (2008) and in high-mass YSOs by Gerakines 
et al. (1999), we conclude that all of our outflow YSOs are 
at an earlier evolutionary stage such that the ice is still cold. 
The cold ice may be the reason why we do not see gaseous 
C2H2, HCN, and CO2, as mentioned above. 

In Fig. 21(a) we plot the optical depth of the 15.3 urn 
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Table 4. Maximum silicate and ice optical depths. For each YSO in an outflow source, this table 
gives the optical depth of the silicate feature at 9.6 um, tq.q, the optical depth and error of the 
polar component of CO2 ice, T p _co 2 > a * ~ 15-3 nm, the optical depth and error of the sum of the 
apolar components of CO2 ice, T ap _co,i at ~ 15.1 [im (see text), the ratio of the apolar to the 
polar components, ap/p, and the normalized central depths and errors measured by simultaneous 
Gaussian fits to the 6.0 and 6.8 (im absorption features, Tg.o and Tg.g, respectively (see text). 
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Figure 21. Optical depths of the polar CO2 component of the 
CO2 ice feature plotted vs. the 9.6 urn silicate optical depth. The 
sources are identified by symbol and by colour as given in the 
figure. Panel (a): a linear relation between the optical depths of 
the two compounds is given by the straight lines, computed by 
least-squares. Panel (b): for the two sources with intercepts clearly 
not at the origin, the excess silicate optical depth is subtracted: 
Ar 9 . 6 = 2.03 for G333.184 and Ar 9 . 6 = 2.08 for G333.466. This 
excess silicate optical depth can be attributed to the intervening 
diffuse ISM, whose dust grains lack CO2 ice (e.g., Whittet et al. 
2007). 

CO2 ice feature (the strong polar component) versus the op- 
tical depth T9.6 of the 9.6 [im silicate absorption feature as 
estimated by our pahfit computations (Section 3.1). The 
plot shows that there is a correlation for individual sources, 
although there are differences between sources, as shown 
by the fits to the CO2 optical depths versus silicate op- 
tical depth. In Fig. 21(b) we subtract the silicate optical 
depth corresponding to the intersection of the fitted line 
with zero CO2 optical depth (rg.e = 2.03 for G333.184 and 
rg.6 = 2.08 for G333.466) - that is, we assume that the 
differences between the sources is due to the foreground ex- 
tinction through the diffuse ISM. Now we see that there is 
much better agreement between sources. We conclude that 



in massive YSOs with ices, the CO2 ice absorption is pro- 
portional to the silicate absorption once the foreground ISM 
absorption, which lacks CO2 ice, is removed. 

3.3.2 The 15 \i.m CO2 ice feature compared to the other 
ice features 

Because of the PAH emission present in every spectrum, the 
features at 6.0 and 6.8 u.m cannot be analysed with the same 
precision as the 15.2 urn CO2 ice feature. In fact, we are not 
able to do any more than determine the normalized central 
depths, T, of the stronger features, which we estimated by 
fitting Gaussians to the two features. We do not fit templates 
to the 6.0 and 6.8 urn features because both are thought to 
be composites of at least two molecules each (e.g., Boogert et 
al. 2008). The measured central depths are given in Table 4 
for the outflow YSOs. 

From the table there is clearly a difference, source to 
source, between each of the ices: CO2, the 6.0 urn feature, 
and the 6.8 u.m feature. The best correlation of the strength 
of the 6.0 u.m feature relative to the CO2 optical depth is 
with the wavelength of the peak of the SED (Table 3 and 
Fig. 19): colder SEDs have a larger 6.0 urn feature relative to 
the CO2 ice optical depth. There is no apparent correlation 
with source morphology or luminosity. 

3.4 Sources with forbidden line emission 

As was mentioned above, there are forbidden emission lines 
present in almost all the spectra. In many of the spectra, for 
example in the background spectra, these forbidden lines 
clearly arise from the diffuse ionized ISM in the G333 cloud. 
On the other hand, some of the outflow YSOs have min- 
ima of the forbidden lines at the YSO, probably because of 
their extreme extinction (Table 3). Consequently, we iden- 
tify sources producing ionizing photons only when there is 
a maximum of the forbidden lines uncorrected for extinc- 
tion at the source location. The sources producing ionizing 
photons are listed in Table 5 and contours showing the con- 
centration of forbidden lines are plotted in Figs. 22 - 26. 

There are two outflow sources (G333.131 and G333.466) 
in Table 5 whose contours show the local production of 
forbidden lines. However, in neither of the outflow sources 
are the contours of maximum forbidden line emission lo- 
cated near the outflow source but are located at some dis- 
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Table 5. Sources producing ionizing photons. The spectrum of each source with localized forbidden lines after background sub- 
traction was analysed as an Hn region. This table gives the electron density, N e , estimated from the ratio of the [S III] 18.7 and 
33.5 urn lines; ratios of neon, sulphur, and hydrogen ions, the ratio of the total neon/sulphur abundance; the radio flux at 5 GHz 
predicted from the strength of the H 7-6 recombination line, £5; and the numbers of hydrogen- ionizing photons, Ni, yc , estimated 
from the neon forbidden line fluxes (see text). Exponents are written as follows: 1.8e48 equals 1.8 X 10 48 . 
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Figure 22. Images of G332. 800— 0.595 with emission line contours. Panel (a): the [Nell] 12.8 urn line contours are plotted in white and 
the [Nem] 15.6 urn line contours are plotted in black on the IRAC band 3 (5.8 um) image. Panel (b): the [S III] 18.7 urn line contours 
are plotted in white and the [S III] 33.5 um line contours are plotted in black on the IRAC band 4 (8.0 urn) image. Panel (c): the [Aril] 
7.0 |im line contours are plotted in white and the [Sin] 34.8 um line contours are plotted in black on the MSX 21 um image. Panel (d): 
the H2 S(2) 12.3 urn line contours are plotted in white and the H2 S(l) 17.0 um line contours are plotted in black on the MIPS 70 um 
image. 



Massive young stellar objects in G333. 2—0.4 17 



2 20 



o 

O 
LO 



40 



o 
o 

8 60 



I 80 



G333.1 63-0.1 00 





_i , , , i_ 



_l , , , L 



_l , , , I , , , I , , , L 



46 44 42 40 38 46 44 42 40 38 
RA (J2000) 1 6 h 1 9 m RA (J2000) 1 6 h 1 9 m 

Figure 23. Images of G333. 163— 0.100 with emission line contours. Panel (a): the [Nell] 12.8 |xm line contours are plotted in white and 
the [No in] 15.6 urn line contours arc plotted in black on the IRAC band 3 (5.8 |xm) image. Panel (b): the [Sin] 33.5 urn line contours are 
plotted in white and the [Sill] 34.8 urn line contours are plotted in black on the IRAC band 4 (8.0 |rm) image. The white cross marks 
the location of the 6.7 GHz methanol maser (Caswell 2009). 



tance, at least 20 arcsec but possibly further since the maps 
don't quite cover the region of forbidden line emission. The 
G333.131 Hn region is not imaged in this group of fig- 
ures because the forbidden lines measured in most telescope 
pointings do not have good enough S/N - the only telescope 
pointings with significant forbidden lines are in the bright 
fuzzy region north of Dec -50° 40' 35" (Fig. 2e). These two 
sources both have the appearance of clusters in their 24 (im 
MIPS images (Figs. 2e and 26d) and thus the forbidden line 
emission may be associated with some star other than the 
YSO producing the outflow. In the G333.466 source, the 
forbidden lines have their maxima at the west side of the 
image, but there is a secondary maximum at the northern 
tip of the outflow seen in the 4.5 u.m IRAC image (Fig. 
26c). In Table 5 and Fig. 27 these regions are described as 
G333.466-0.164W and G333.466-0.164N, respectively. 

We analysed the forbidden lines and the H 7-6 recombi- 
nation line using the code described by Simpson et al. (2007) 
with updates from Table 2 for an assumed electron temper- 
ature, T e , of 7000 K. The line fluxes are given in Table A3. 
The results are also in Table 5 for electron density, N e , es- 
timated from the ratio of the [S ill] 18.7/33.5 [im lines, and 
the ratios of various neon, sulphur, and H + ions. The IR 
forbidden lines have very low temperature sensitivity. How- 
ever, the H + recombination line is sensitive to T e such that 
the ionic abundances with respect to H + are approximately 
proportional to T~ . Thus the fact that the Ne/H and S/H 
abundances are generally not very different from what one 
would expect for sources closer to the Galactic Centre than 
the Orion Nebula (Orion Nebula Ne/H = 1.01 x 10 -4 and 
S/H = 7.7 x 10~ 6 , Rubin et al. 2011), given the Galactic 
abundance gradients (e.g., Simpson et al. 1995) and the 
abundances in G333.6-0.2 (Ne/H = 3.4 x 10~ 4 and S/H 
= 1.3 x 10~ 5 , Simpson et al. 2004), leads us to infer that 



T e = 7000 K is a reasonable value. Such an electron temper- 
ature is close to that expected for photoionized H n regions 
with these abundances and much lower than the tempera- 
tures predicted for gas ionized by shocks. 

We also used the observed H 7-6 fluxes to predict the 
radio flux at 5 GHz, S5, and the observed [Nell] 12.8 u.m 
flux to estimate the number of ionizing photons, A?Lyc, us- 
ing equation 1 of Simpson & Rubin (1990). These parame- 
ters can also be estimated from the observed forbidden line 
fluxes, which is useful because the forbidden lines are much 
stronger than the radio flux or IR hydrogen recombination 
lines, both of which are often not detectable or detectable 
only with very large uncertainties. 

We start with equation 1 of Rubin (1968) but multiply 
the left side by [1 + f t < He+/(H+ + He+) >] (as was done 
by Simpson & Rubin 1990), to get 



Lyc 



He 4 



(H+ + He+) 



-i: 



(1) 



that is, the number of stellar photons in the Lyman con- 
tinuum able to ionize the hydrogen in an Hn region, Ni, yc , 
augmented by the ionizing photons from He + recombina- 
tion, is equal to the total number of recombinations to neu- 
tral H and He (ignoring heavy elements) integrated over 
the whole Stromgren sphere of radius Rs and volume V, 
where N e is the electron density and /; (« 0.65) is the frac- 
tion of helium recombination photons to excited states which 
are energetic enough to ionize hydrogen (Simpson Sz Rubin 
1990). For T e from 5000 K to 10,000 K and log N e ~ 2.5, 
the total recombination coefficient for Case B, qb, equals 
~ 2.587 x 10" 13 i" - 81 cm 3 s" 1 for H (Storey & Hummer 
1995) and ~ 2.755 x 10~ 13 t~ ' 78 cm 3 s" 1 for He (Hummer 
& Storey 1998), where t = T e /10 4 . 

We consider some heavy element X, with abundance 
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Figure 24. Images of G333. 340— 0.127 with emission line contours. Panel (a): the [Nell] 12.8 um line contours arc plotted in white and 
the [Ncin] 15.6 um line contours are plotted in black on the IRAC band 3 (5.8 um) image. Contours of the [Siv] 10.5 um line look 
similar to the [Ncm] line contours. Panel (b): the [S III] 18.7 um line contours are plotted in white and the [S III] 33.5 um line contours 
are plotted in black on the IRAC band 4 (8.0 um) image. Contours of the [Fein] 22.9 um line look similar to the [S III] 33.5 um line 
contours. Panel (c): the [Sill] 34.8 um line contours are plotted in white on the MIPS 24 um image. Contours of the [Fell] 26.0 um line 
look similar to the [Sill] line contours. The centre of the MIPS image is saturated. Panel (d): the H2 S(2) 12.3 um line contours are 
plotted in white and the H2 S(l) 17.0 um line contours are plotted in black on the IRAC band 4 (8.0 um) image. 



with respect to hydrogen Ax /Ah, which has an ion- 
ized state with ionic density AT, that emits a collision- 
ally excited forbidden line with volume emissivity e, e — 
hvA{N up per/Ni)Ni, where h is Planck's constant, v is the 
frequency of the transition, A is the Einstein transition prob- 
ability, and N upper is the density of the ions in the upper 
level of the transition. For computational convenience we 
define e such that its only dependencies are N e , T e , and the 
atomic physics parameters (Table 2) . With that we have 



6 Ax A H Ae ' 



(2) 



where A p is the proton density. In our code we divide e by 
because the emissivity at low densities is proportional 
to and the resulting e' /N% has little density dependence 
until the density is high enough to produce collisional de- 
excitation. This collisional de-excitation allows estimation 
of A e from the ratios of two lines from the same element, 
such as S ++ (e.g., Simpson et al. 2004). 



The flux from this line, observed over the whole emitting 
region, is given by 



4ttD 2 J A| A x A H A e e 



(3) 



where D is the distance to the Hn region, which we assume 
is fully ionized (N p — Ah). Assuming uniform excitation 
and substituting the J N^dV derived from eq. (3) into eq. 
(1), we get 



2.587 x 10 _13 ^°- 81 47rD 2 F 



4 



1 + fi ( (H+ H +He+) 



(4) 



where the brackets on < Ni/Nx >, etc., indicate an average 
over the emitting volume (see Simpson et al. 1995, 2004, for 
further examples of the notation) , and the lower limit arises 
because of the possibility that the measured flux originates 
from less than the total ionized volume. Here we have as- 
sumed the as for H since the correction for the different 
value for He as given above is so small. 
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Figure 25. Images of G333. 375— 0.202 with emission line contours. Panel (a): the H2 S(2) 12.3 |xm line contours arc plotted in white 
and the [Nell] 12.8 |xm line contours are plotted in black on the IRAC band 3 (5.8 um) image. Panel (b): the [Sin] 34.8 |xm line contours 
are plotted in white and the [S III] 33.5 um line contours are plotted in black on the IRAC band 4 (8.0 |xm) image. 



For use with forbidden lines (when no reliable hydrogen 
recombination line or radio flux is available), it is prefer- 
able to use a line that gives the least amount of uncer- 
tainty - that is, low extinction, ionization correction fac- 
tor (< Nx/Ni >) close to unity, and relatively well-known 
abundance (< Nx/N K >). For the MIR observations of low- 
excitation nebulae (no He + ) discussed in this paper, a good 
line is the 12.8 [im line of Ne + , which has very little den- 
sity sensitivity. For the line parameters and references in 
Table 2 we calculate e'/A" 2 w 9.38 x lO" 2 ^- 276 erg cm 3 
s" 1 for log N e = 2.5 (N e in cm" 3 ) and T e between 5000 K 
and 10000 K. This gives 



A^Lyc ^ 3.30xl0 55 £ 



-0.525 n 2 



kpc-^12.8 



Ne + 



/Ah 

\ Anc 



(5) 



where F12.8 is the observed flux of the [Nell] 12.8 u.m line 
in W m -2 and Dk pc is the distance in kpc. 

In Table 5 we estimate A^yc from the observed neon 
fluxes and predict radio fluxes from the H 7-6 recombina- 
tion line flux. Here we assumed a Ne/H ratio of 1.5 x 10~ 4 
(the weighted average of the Ne/H abundance ratios in Ta- 
ble 5), T e = 7000 K, and N e = 316 cm~ 3 , and used the 
volume cmissivities computed by the forbidden line analysis 
program. Since there are ionized lines in all spectra from the 
background, background fluxes were first subtracted. These 
background fluxes may be overestimates since they were 
taken from positions relatively close to the sources where 
the intensity is at a minimum (see the plots of the SH slit 
locations in Figs. 8 - 18), and consequently the estimates 
of S5 and N-Ly C could be underestimates. This is particu- 
larly true of the two sources with the smallest AL yc , where 
the background intensity is a sizable fraction of the average 
source intensity. 

Assuming, then, that the sources can be analysed as 
photoionized H n regions, we compare the estimated ionic ra- 



tios to ratios that we compute with CLOUDY, version 08 (Fer- 
land et al. 1998). We computed models assuming that the 
ionizing sources can be described as zero-age main-sequence 
(ZAMS) stars, since such stars have the hottest effective 
temperatures, T e s for their luminosities. We took the stel- 
lar radius-T c ff -luminosity relation for ZAMS stars from the 
'canonical' Z = 0.02 pre-main-sequence evolution tracks of 
Bernasconi & Maeder (1996) and used both black bodies 
and the log g — 4 atmospheres of Sternberg et al. (2003) for 
T cff = 25000 K to 36000 K or 40000 K. cloudy was used to 
compute A^Lyc vs. T c ff for each atmosphere. 

In Fig. 27(a) we plot the computed A^yc vs. luminosity 
for the ZAMS stars along with the results for our sources 
from Table 5. Luminosity is plotted as the ordinate and 
Al yc as the abscissa because A^yc can be thought as a proxy 
for T e ff. Upper limits for the non-ionizing sources were es- 
timated from the neon flux uncertainties. We also plot the 
N-Lyc vs. luminosity relation for luminosity class V stars from 
Vacca et al. (1996), Sternberg et al. (2003), and Martins et 
al. (2005) in Fig. 27(a). Because our sources lie in the region 
of the ZAMS stars on the plots, we again conclude that it 
is likely that the sources are ionized by stellar photons and 
not shocks. 

In Fig. 27(b) we plot the observed Ne ++ /Ne + ratios vs. 
the FIR luminosities and the ratios predicted by the same 
models. It is notoriously difficult to get models that pre- 
dict the observed Ne ++ /Ne + ratios because of the extreme 
sensitivity of the Ne ++ ionization level (which requires 41 
eV) to details of the input stellar atmosphere models (e.g., 
Simpson et al. 2004; Rubin et al. 2008). In general, our ob- 
servations indicate higher Ne ++ /Ne + ratios than are pre- 
dicted by the ZAMS star models (Hn region models using 
luminosity class V star atmospheres would have even lower 
Ne ++ /Ne + ratios at given luminosities because of the lower 
stellar T e ff). Blackbody atmospheres were included for this 
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Figure 26. Images of G333. 466— 0.164 with emission line contours. Panel (a): the H2 S(2) 12.3 um line contours are plotted in white 
and the H2 S(l) 17.0 urn line contours are plotted in black on the IRAC band 2 (4.5 um) image. The white cross marks the location of 
the 6.7 GHz methanol maser (Caswell 2009); this is also the location of the YSO with the ice absorption features. Panel (b): the [Sill] 
34.8 um line contours arc plotted in white and the [Fen] 26.0 urn line contours are plotted in black on the IRAC band 4 (8.0 um) image. 
Panel (c): the [Nell] 12.8 um line contours arc plotted in white and the [Nein] 15.6 urn line contours are plotted in black on the IRAC 
band 2 (4.5 um) image. The [Aril] 7.0 um line contours look similar to the [Nell] contours. Note the displacement of the ionized lines 
from the bright outflow emission, but there is a secondary contour maximum off the northern tip of the outflow. Panel (d): the [S in] 
33.5 um line contours are plotted in white and the [S III] 18.7 um line contours are plotted in black on the MIPS 24 urn image. Note that 
the [S in] line contours peak at the brightest (saturated) pixels of the MIPS image. 



project because their use in Hn region models gives better 
agreement with observations than other model atmospheres. 
Most of the sources plotted in Fig. 27(b) have Ne ++ /Ne+ 
ratios lower than predicted by the blackbody models; how- 
ever, the fact that two sources, G333.340 and the Hn region 
north of G333.131, have higher Ne ++ /Ne + ratios than even 
those predicted by the blackbody models may indicate an 
additional source of ionization, such as shocks. This will be 
discussed in Section 3.6. 



3.5 Molecular hydrogen lines 

It is not simple to derive the PDR or shock properties from 
the observed lines for several reasons: (1) The main indica- 
tors in a general sense are the ratios of the H2 rotational 



lines, but only the first three [S(0), S(l), and S(2)] are read- 
ily detectable with the SH and LH modules' \/SX = 600 
resolution, since the SL module's resolution is only 64 — 128. 
(2) The H 2 S(3) line falls in the bottom of the 9.6 um sil- 
icate absorption feature and the S(5) line is blended with 
the usually strong [Aril] 6.98 (im line. (3) Other lines, 
such as [Fell] 26.0 (im and [Sill] 34.8 um, which are im- 
portant PDR indicators (Kaufman, Wolfire, & Hollenbach 
2006), are also found in Hn regions. Moreover, one must 
assume some depletion factor onto grains. (4) Finally, and 
most importantly, all the above mentioned lines also arise 
in shocks. Considering the first two items, it is perhaps not 
surprising that we detect the S(3), S(5), and S(7) lines in 
only one source, G332. 813— 0.700, and there very poorly be- 
cause of the strong continuum (see Fig. 3b), and no S(4) or 
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Figure 27. Panel (a): source luminosity plotted vs. the numbers of Lyman continuum photons per second for the six sources producing 
their own ionizing photons and 3 cr upper limits for the nine YSOs without measurable H II regions. Estimated luminosities and Ni^y C 
for O stars of luminosity class V as determined by Vacca, Garmany, & Schull (1996), Sternberg, Hoffmann, & Pauldrach (2003), and 
Martins, Schaerer, & Hillicr (2005) are also plotted. The connected diamonds and plus signs are models for Hn regions ionized by ZAMS 
stars (lower luminosities than luminosity class V stars) with T e g ranging from 25000 K to 36000 K for black body SEDs and T e g ranging 
from 32000 K to 41000 K for SEDs from Sternberg et al. (2003). The source luminosities are the average of the best fitting models 
from the SED fitter and the uncertainties are the standard deviations of these fits. An additional uncertainty due to uncertainty in the 
distance ~ 10% is not plotted. Panel (b): source luminosity plotted vs. the neon ionic ratio for the seven Hn regions (note that there 
are two positions for the G333.466 Hll region, see text). 



S(6) at all. We also detect the S(5) and S(7) lines in the 
G333.466-0.164 outflow, but no S(3). 

The usual way of distinguishing shocks from PDRs is 
that shocks do not have much thermal continuum (e.g., Hol- 
lenbach, Chernoff, & McKee 1989), unlike PDRs where the 
grains are heated by the stellar photons. We cannot use this 
criterion here because of the strong FIR continuum from the 
YSOs. Instead, we look at the morphology of the line maps, 
as we did for the ionized lines, and also compare the Hb line 
ratios to theoretical models of both shocks and PDRs. Con- 
tour plots of PDR/shock lines in outflow sources are shown 
in Figs. 28 - 30. We see no correlation of any of these lines 
with the outflows. We conclude that the morphology seen 
in the images is most likely due to the variations in extinc- 
tion that are present in these sources. The best agreement is 
with the morphology of the PAH emission seen in the IRAC 
8.0 |im images, as would be expected if the H2 lines are 
formed in PDRs like the PAHs. 

Compared to the models, we typically see that the 
sources have more H2 S(0) relative to S(l) than the shock 
models of Hollenbach & McKee (1989) and Flower & Pineau 
des Forets (2010) and more H2 S(2) relative to S(l) than the 
PDR models of Kaufman et al. (1999, 2006) 5 . We consider 
only the shock models with velocities less than 40 km s _1 for 



5 http://dustem.astro.umd.edu/index.html 



the outflow YSOs with localized H2 because at higher veloc- 
ities, the models predict detectable lines of ionized elements, 
which are not observed (see the next two sections for possi- 
ble exceptions). The higher H2 S(0)/S(1) ratio suggests that 
there is more lower temperature gas present (at ~ 100 K) 
than produced by shocks, which are generally much hotter. 



The PDR models of Kaufman et al. (1999, 2006) all have 
S(2)/S(l)< 1 (the observed ratio is typically ~ 1) except for 
the extreme model range of very high density, n (n ~ 10 4 ' 5 
cm -3 ), and either very high or very low Go (where Go is 
defined as the incident ultraviolet flux from 6.0 to 13.6 eV 
in units of 1.6 x 10 -3 erg cm -2 s _1 ). The models with the 
S(1)/S(0) ratio in the range of ~ 0.2 to ~ 3 (the majority 
of the data) have either n ~ 10 3 cm~ 3 and Go > 10 2 ' 5 
or n > 10 3 cm~ 3 and Go < 10 2 . We conclude that the 
observed H2 ratios are in reasonable agreement with the high 
density, low Go models of Kaufman et al. (2006), especially 
if we consider multiple components along the line of sight. 
In fact, we should expect to sec PDR H2 line emission from 
our spectra, given the large amount of PAH emission seen 
in our spectra and in the IRAC 8.0 urn images. This PDR 
emission probably originates in the cloud as a whole and 
is not local to the outflow YSOs, which are all seen to be 
embedded in local dust lanes of higher extinction. 
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Figure 28. Image of G332. 813— 0.700 with emission line contours. Panel (a): the H2 S(l) line contours are plotted in white and the H2 
S(0) line contours are plotted in black on the IRAC band 2 (4.5 urn) image. Panel (b): the [Fen] 26.0 um line contours are plotted in 
black on the IRAC band 4 (8.0 um) image. 
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Figure 29. Image of G332. 942— 0.686 with emission line con- 
tours. The [Sill] 34.8 um line contours are plotted in white and 
the H2 S(l) line contours are plotted in black on the IRAC band 
2 (4.5 um) image. It is curious that these contours appear or- 
thogonal, because it usually expected that both lines are formed 
in PDRs, that is, in the same gas. The white cross marks the 
location of the 6.7 GHz methanol maser (Caswell 2009). 



3.6 Shocks 

Both the red sources and the outflow sources show line emis- 
sion indicative of shocks. Although the peak of the low exci- 
tation forbidden lines is at the location of the YSO in the red 
sources G332.800 and G333.340, in both of these sources the 
high excitation forbidden lines peak well off the YSO (Figs. 




Figure 30. Image of G333. 184—0.091 with emission line con- 
tours. The [Sill] 34.8 urn line contours are plotted in white and 
the H2 S(0) line contours are plotted in black on the IRAC band 
4 (8.0 um) image. The white cross marks the location of the 6.7 
GHz methanol maser (Caswell 2009). 



22 and 24). This is the opposite of what is expected for 
an H 11 region, where ions requiring higher energy photons 
for photoionization have much smaller Stromgren spheres 
immediately surrounding the ionizing star than ions need- 
ing the much more abundant photons with energies only a 
little larger than 13.6 eV. This is particularly apparent in 
the [Neil] 12.8 |im and the [Nem] 15.6 |J.m linos (panel a 
of both figures), which have similar extinction and critical 
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densities. Ions of intermediate ionization potentials, such as 
[S in] (panels b of Figs. 22 and 24), have intermediate mor- 
phologies. We suggest that the offset high excitation lines 
are produced in an outflow that is shocking the surrounding 
envelope or ISM. The alternative explanation, that a wind 
has blown out a cavity and the stellar photons are ionizing 
the cavity edges, is not in agreement with the observation 
that the highest excited lines ([Nem] and [Siv]) are found 
farther from the star than lower excitation lines like [Sin]. 

The various line ratios can be used with models to esti- 
mate the shock velocities, given the other necessary model 
parameters such as density and magnetic field, neither of 
which we know. For these two red sources, the most obvious 
shocked lines are [Nem] 15.6 [im and [Siv] 10.5 urn (sim- 
ilar morphology to the [Nem] lines); the lower excitation 
lines like [Nell] 12.8 |im, [Aril] 7.0 |im, [Fell] 26.0 urn, and 
[Sill] 34.8 u.m have extensions along the line of the [Nem] 
emission. The upper limit to the velocity is strongly con- 
strained by the fact that there is no significant emission of 
the higher excitation lines of [O iv] 25.8 urn or [Ne V] 14.3 
or 24.3 um. Only H2 S(2) appears to be affected by the 
shock in G333.340 (Fig. 24d) but no H2 lines are affected in 
G332.800 (Fig. 22d). 

We find in the literature J-shock models for shock ve- 
locities of 30 - 150 km s _1 (Hollenbach & McKce 1989) and 
100-1000 km s _1 (Mappings III, Allen et al 2008). None of 
these models agrees with the data for G333.340 because all 
shock models with the [Nem] 15.6 um line stronger than the 
[Nell] 12.8 urn line also have strong enough [Oiv] 25.85 urn 
emission that we would have detected the line. Consequently, 
we must estimate shock velocities near the bottom of the 
range of model velocities, that is, ~ 100 km s _1 , but em- 
phasize that this value has a large uncertainty because we do 
not know the extremely important magnetic field parameter 
and because of the lack of fit. Further modelling is beyond 
the scope of this paper. 

G333.163 and G333.375 also have the emission extended 
away from the YSO, although it is much less definite that 
the extension is due to an outflow (Figs. 23 and 25). 

The shock velocities for the outflow sources are much 
lower. These sources have no forbidden line emission from 
any ions requiring as much as 13.6 eV for ionization (al- 
though perhaps there is some emission beyond the end of 
the outflow in G333.466, Fig. 26c). For these sources the 
outflows, as identified by the excess IRAC 4.5 urn emission, 
are where we detect the [S 1] 25.25 urn line. The [S 1] 25.25 urn 
line is found only in shocks as the result of shock dissociation 
of molecules - the ionization potential of S + is only 10.4 eV, 
such that any atomic sulphur is always ionized in the diffuse 
ISM. The shock models of Hollenbach & McKee (1989) pro- 
duce substantial [S 1] 25.25 u.m line emission, which is seen 
in both YSO outflows (e.g., Haas, Hollenbach, & Erickson 
1991) and in supernova remnants (e.g., Neufeld et al. 2009). 

Estimates of the maximum surface brightnesses of the 
five sources with detected [S 1] 25.25 u.m lines are given in 
Table 6. This line can be difficult to detect because of the 
small line/continuum ratio at the Spitzer resolution and the 
substantial (< 5 per cent) fringing at this wavelength in the 
LH spectra with strong continua. For two of the sources, the 
line was detected only in the central position and only after 
the spectra were extracted and defringed with smart. For- 
tunately, the [S 1] 25.25 urn line occurs in two separate orders 



Table 6. Maximum observed intensities of the [Si] 
25.25 urn line. The [Si] lines in G332. 725-0. 620 and 
G332. 813— 0.700 were measured with SMART for the 
centre positions in the LH map (LH aperture is 247.5- 
arcsec) because they needed defringing. The [S 1] lines 
in the other sources were extracted with CUBISM with 
smaller extraction apertures in order to estimate the 
maximum intensities, which occur in the outflows. 



Source 




Position 


Intensity 










(W m -2 sr~ 


- 1 ) 


G332.725- 


-0.620 


Centre 


8±1 x 10" 


9 


G332.813- 


-0.700 


Centre 


4±1 x lO- 


9 


G332.942- 


-0.686 


Outflow 


ll ±2 x 10" 


-9 


G333.131- 


-0.560 


Outflow 


7±1 x 10" 


9 


G333.466- 


-0.164 


Outflow 


18 ±4 x 10" 


-9 



of the LH array, and thus we can require that we measure 
the line in both orders to deem that we have a detection. 
For the other three sources, the outflow regions have much 
less continuum than the YSO positions; consequently, fring- 
ing is less important and the lines are easily seen in spectra 
extracted with CUBISM. 

We do not see any increase in H2 emission at the lo- 
cation of the outflows and [S 1] emission, although generally 
the H2 line ratios agree with neither the PDR nor the shock 
models. If we assume that the H2 S(0) line arises in PDRs 
and the S(l) and S(2) in shocks, we can estimate that the 
shock velocity is between ~ 10 km s _1 and ~ 40 km s" 1 
for the S(2)/S(l) line ratio ~ 1, using the models of Flower 
& Pineau des Forets (2010). The velocities cannot be higher 
than this because the models of Hollenbach &: McKee (1989) 
produce lines of [Nell] 12.8 [im, which we do not detect in 
the outflow sources. 

When the excess IRAC 4.5 [im emission was first dis- 
covered, an immediate suggestion was that it could be due to 
high-J H2 emission from shocks in the outflows (e.g., Smith 
et al. 2006; Cyganowski et al. 2008; Chambers et al. 2009). 
In fact, De Buizer & Vacca (2010) found that this must be 
the case for G19.88— 0.53, which has very little continuum 
at 4.5 um and substantial H2 emission. (Their other source, 
G49.27— 0.34, has only continuum at 4.5 [im and no H2 emis- 
sion.) However, this cannot be the case for our sources. In 
Table 7 we list the continuum intensities that we have mea- 
sured for the excess 4.5 urn emission positions in the IRAC 
band 2 images and measurements or upper limits for the H2 
S(7) line at 5.5 urn. In Table 8 we list our measurements of 
the H2 lines that we measured in G332.813, the only source 
with detectable S(3) - S(7) lines. Even though we have mea- 
sured only the H2 S(l) to S(7) lines and the lines in the 
IRAC band 2 bandpass are S(9) and S(10), the measured 
4.5 yim continuum is too large by a factor of at least 10 to 
be due to H2 line emission. We conclude that some other 
source of emission must be the dominant contributor to the 
IRAC band 2 intensities. This could be either CO emission 
(e.g., Marston et al. 2004) or perhaps scattered light from 
the large grains that are sometimes found in the cores of 
star- forming regions (Pagani et al. 2010). Some scattered 
light is almost certainly present, since the excess 4.5 nm 
emission regions are clearly present in the IRAC band 1 im- 
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Table 7. IRAC surface brightnesses and H2 intensities. In response to the suggestion that the IRAC Band 2 extended 
emission could be due to emission from H2 in outflows, this table lists the IRAC Band 2 intensities in both outflow regions 
and nearby background positions, and the difference that is attributable to the outflow in units of both MJy sr -1 and 
W m~ 2 sr -1 . The IRAC bandwidth for Band 2 is 1.015 [ira at effective wavelength 4.493 urn. The last column lists the 
measured H2 S(7) intensities or 3<r upper limits at the same positions. Since the background was not subtracted, the H2 
S(7) intensities arc upper limits to the H2 S(7) intensities in the outflow region, especially since there is no or little apparent 
increase in the H2 intensity at the outflow positions in the table. The H2 intensites at 4.5 mil [H2 S(9) and S(10)] could be 
slightly but not substantially more than the H2 S(7) intensity, given the flat extinction law at these wavelengths (Indebetouw 
ct al. 2005). Exponents are written as follows: 2.43e-6 equals 2.43 X 10~ 6 . 
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ages at 3.6 ^irn and also in the 2MASS Ks and sometimes 
H images. 

3.7 Individual sources 

In this section we comment on individual sources. 

(i) G332. 725-0. 620. This source is missing the SL first 
order spectrum (SL1) from 7.4 - 14 \irci because the pre- 
dicted saturation of the IRS Peak-Up arrays, which share 
the same detector array as the SL spectrograph module, 
would render all spectral pixels in the same rows on the 
detector array uncalibratable. This is unfortunate because 
the SL1 wavelength range provides the best estimate of the 
depth of the 9.6 u.m silicate feature. Hence our estimates of 
the extinction in this source are more uncertain than the 
estimates for the other sources. 

(ii) G332. 800-0. 595. This source is too extended to 



map with the IRS slits touching. The result is that there 
probably is missed structure in the maps and the position 
of the ionized line outflow is poorly defined (see Fig. 22) . It 
is also not known which of the two bright sources seen in the 
centre of the IRAC images (see Fig. 22 and Fig. 31 of the 
Supporting Information) is the YSO. The source on the left 
is extended and is bright in 2MASS J, H, and Ks bands. 
The source on the right is compact and much redder, being 
only barely visible in the 2MASS J band. The sources are 
approximately 10 arcsec apart, which means they cannot be 
distinguished in the 19 arcsec resolution MSX 21 u.m image 
(Price et al. 2001) of Fig. 22(c) (the MIPS 24 urn image was 
saturated) or the MIPS 70 u.m image. The MSX Catalog 
gives the coordinates for the source as RA 16 h 20 m 16!6, 
Dec —50° 56' 23", that is, the left-hand source. 

(iii) G332. 813- 0. 700. This source is the only outflow 
source without a methanol maser but it is also the only 
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Table 8. H 2 lines in G332.813-0.700. The H 2 lines were well-detected only in G332.813-0.700. This 
table gives the line intensities at the YSO centre and at the position of the strongest lines, 'north' 
(the north position is at RA 16 h 20 m 47!6, Dec. -51°00'6"and the centre position is at RA 16 h 20 m 
48? 1, Dec. — 51°00'15"). The line intensities were corrected for extinction (/ cor r) at both positions 
with the measured T9.6, but the relative intensities for the centre position (/ C orr centre) show that 
the correction is much too large for the line with the largest correction, S(3). The last column gives 
the corrected intensities for a much lower value of Tg.g; we infer from the more reasonable line ratios 
that the H2 lines originate from a region foreground of the YSO in the Center. PDR parameters, n 
(density in cm -3 ), and Go (incident ultraviolet flux from 6.0 to 13.6 eV in units of 1.6 X 10 — 3 erg 
cm~ 2 s _1 ), were estimated from the corrected line ratios for these two positions (models of Kaufman 
et al. 2006). Exponents are written as follows: 1.88c-8 equals 1.88 X 10~ s . 
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source with detectable H 2 S(3) lines (Table 8). The H 2 lines 
peak in the PDR (region of strong PAH emission) to the 
north-west of the YSO and not in the region of extended 
IRAC band 2 (4.5 pm) emission as seen in Fig. 28(a). The 
[Fell] 26 pm line also peaks in the PDR regions that show 
strong PAH emission at 8 pm, as seen in Fig. 28(b). We 
suspect that this PDR region is somewhat foreground of 
the main source (which lies in northeast-southwest lane of 
higher extinction) because the extinction-corrected S(3) flux 
is much higher than the fluxes of the S(l) and S(5) lines (the 
S(3) lines is at 9.66 pm). Other YSO spectra may appear 
to have the S(3) 9.66 pm line, but there are numerous bad 
pixels at this wavelength in the SL array; since a resolution 
element in the IRS is 4 pixels and since the H 2 lines should 
be extended and not point-like, we do not believe we have 
identified a line unless we can see it on the array in a number 
of contiguous pixels, as we can for this source but not the 
others. 

(iv) G332.942~0.686. In this outflow YSO the H 2 line 
emission peaks to the northwest, in the direction of the PAH 
emission as seen in the IRAC 8.0 pm image, but the [Sill] 
34.8 um line emission increases to the southwest. Perhaps in 
this case the Si + is found in the diffuse Galactic ISM but the 
H 2 is more associated with the G333 GMC, located towards 
the northwest. These contours are plotted in Fig. 29. 

(v) G332.963-0.679. This outflow YSO is the most 
luminous source of the little cluster that includes the red 
sources G332. 967-0.683 and G332. 960-0. 681. We see in 
Supporting Information Fig. 32(e) that it is the strongest 
source in the MIPS 24 pm image, and in fact, all the longer 
wavelength images (MIPS 70 pm, and Herschel PACS and 
SPIRE 150-600 urn images) peak at its location. The ex- 
tinction map (Fig. 12) shows that much of the extinction in 
the region is due to the extended envelope (~ 20" « 7 x 10 4 
au). The red sources are foreground, as seen by their lower 
extinctions. 



(vi) G332.967- 0.683. This YSO and the similar red 
YSO G332. 960— 0.681 cannot be detected in images at wave- 
lengths longer than the MIPS 24 pm image, where they are 
both much fainter than the outflow YSO, G332. 963-0.679. 
At least part of the reason for lack of detection is spatial res- 
olution for the Spitzer 70 pm image and the Herschel long 
wavelength PACS and SPIRE images. As a result of lack of 
measurement at the longer wavelengths, the SED and total 
luminosity are poorly measured (Fig. 19). 

(vii) G333. 125-0. 562. This is the name Garay et al. 
(2004) gave to the first detection of this massive, dense 
cold core. We included it in our survey since its SiO 
emission (Lo et al. 2007) is thought to arise from shocks 
associated with an outflow. The YSO with the ice ab- 
sorption features, G333. 131-0.560, is GLIMPSE source 
SSTGLMC G333. 1309-00.5601. The other stars seen by 
the IRAC that are also bright at 24 pm are SSTGLMC 
G333.1310-00.5658 and SSTGLMC G333.1295-00.5683. 
Unfortunately, they lie beyond the range of our SL map 
(Fig. 2c), but since they have similar IRAC colours to SST- 
GLMC G333.1309-00.5601, we think it likely that they also 
have ice-coated grains in their surrounding envelopes. The 
source G333. 128—0.560 is probably the 70 pm source since 
our spectral maps show that it is the brightest object at the 
longest wavelengths, 36 pm, even though it is less bright 
than G333.131 at 24 pm. 

The shocked [S 1] 25.25 pm emission is seen in LH tele- 
scope pointings 70 and 71; these slit positions are marked 
on Fig. 2(e) and correspond approximately to the extended 
IRAC 4.5 pm emission (seen more easily in the 3-colour 
image of Cyganowski et al. 2008). It is not possible to iden- 
tify which of the YSOs is the source of this outflow be- 
cause it is faint and because both G333.131 and G333.128 
lie along the same line, as do the 6.7 GHz Class II methanol 
masers (Caswell 2009; Fig. 2b). The SiO emission occurs at 
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G333. 131— 0.560 whereas the 4.5 urn and [Si] outflow region 
coincides with emission of HCO+ (Lo et al. 2011). 

We did not have enough flux data to model both 
G333.131 and G333.128 separately with the online SED fit- 
ter of Robitaille et al. (2007) - the results in Table 3 are 
from a composite of both sources. A model with a larger 
inclination angle, larger model age, and smaller envelope 
mass would have properties similar to G333. 128— 0.560 - 
more flux at 35 - 70 u.m but much less flux at 24 urn (Fig. 
2e) and 3 mm (Lo et al. 2011). 

The fuzzy blob to the north of G333.131 appears to be 
an H n region, with forbidden lines and substantial 24 [im 
continuum (Fig. 2e). However, it is unusual for an Hn region 
in that there is no obvious continuum emission at the longer 
wavelength maps (Spitzer MIPS and Herschet). If it is ac- 
tually associated with the two stars at its location, 2MASS 
16213573-5040306 and 2MASS 16213662-5040334, it may be 
a foreground object since these stars are not very red and 
the former is seen on the red Digital Sky Survey image. The 
SL spectra show only continuum at 5.2 urn and background 
PAHs. 

(viii) G333. 163-0. 100. The spectral maps of this source 
show extensions to the north-east in the lines mapped by LH: 
[S III] 33.5 urn and [Sill] 34.8 urn, plotted in Fig. 23(b), and 
also [Fen] 26.0 \im and [Fein] 22.9 |xm [this region was not 
covered in the contours of Fig. 23(a) taken from the much 
smaller SH map]. The SL map shows a similar morphology 
for the [Aril] 7.0 urn line. This is somewhat orthogonal to 
the apparent direction of the red extended emission. This is 
the only red source with a 6.7 GHz methanol maser (Caswell 
2009; Fig. 23b). 

(ix) G333. 184- 0. 091. This outflow YSO has a distinct 
absorption lane visible in the IRAC images, especially the 
8 urn image, that is perpendicular to the apparent outflow 
direction (compare Figs. 14 and 30). This morphology is 
still apparent in the 24 um MIPS image (Fig. 33e of the 
Supporting Information). It is much too big to be a disc, 
but could be indicative of a toroidal structure in the envelope 
with a much larger optical depth. The morphology of the H2 
S(0) and [Sill] 34.8 um lines (contours plotted in Fig. 30) is 
probably due to absorption of background gas by the YSO 
envelope. 

(x) G333. 212-0. 105. This red YSO does not appear to 
be present in the MIPS 70 urn image (Supporting Informa- 
tion Fig. 34f), but that is due to the low Spitzer resolution 
and the high background from the other sources in its imme- 
diate vicinity. The source is definitely visible in the higher 
resolution Herschel PACS and SPIRE images, even at the 
longest SPIRE wavelength of 600 |im. 

(xi) G333.340- 0.127. This red source has a distinct 
outflow to the south-west as seen in the higher-excitation 
forbidden lines, particularly [Nem] 15.6 urn and [Siv] 
10.5 urn (Fig. 24). These high excitation lines peak well 
away from the YSO central region, which does include the 
peak of the low excitation lines, like [Nell] 12.8 um, [Sill] 
34.8 urn, and [Fell] 26.0 urn. Lines of intermediate excita- 
tion, like [S III] 18.7 and 33.5 um and [Fein] 22.9 um, have 
intermediate morphologies. The iron lines are not plotted in 
these figures because they are much weaker and noisier, but 
the morphology is clear. There is a clump of warm dust in 
the outflow direction, but closer to the YSO, as seen in the 
MIPS 24 um image and also in the H2 line contours. 



(xii) G333. 375- 0. 202. This red source is too low exci- 
tation to test for outflow morphology from high excitation 
lines. The [Nell] 12.8 urn lines are slightly elongated, similar 
to the structure seen in the PAH emission at 5.6 and 8.0 urn 
and to the H 2 S(0) line (Fig. 25). 

(xiii) G333. 466— 0.164- There are at least two sources 
in this apparent cluster: (1) the outflow YSO with the ice 
absorption features that is the peak at 70 urn and longer 
wavelengths, and (2) an H 11 region located to the west of the 
YSO that is the continuum source seen in the MIPS 24 urn 
image (Fig. 26). The forbidden lines have their maxima at 
the west side of the image, but there is a secondary maxi- 
mum at the northern tip of the outflow seen in the 4.5 urn 
IRAC image (Fig. 26c). In Table 5 and Fig. 27 these regions 
are G333.466-0.164W and G333.466-0.164N, respectively. 
There is also some ionized gas in the outflow (Walsh et al. 
1998). The YSO with the ice features, G333.466-0.164, is 
the location of the 6.7 GHz methanol maser (Caswell 2009) 
and is at the south end of this outflow (Fig. 26a). Fig. 26(a) 
shows that the H2 lines follow the morphology of the PAH 
emission with a big extinction lane seen in the IRAC images. 



4 DISCUSSION 

We review what we have learned about the YSOs in our 
survey of the G333 GMC. Even though the FIR SEDs (Fig. 
19) and luminosities (Table 4) are similar, the spectra of 
the outflow sources are as different from the spectra of the 
red sources as their appearances are different in the IRAC 
bands: (1) The outflow sources have deep ice features at 6.0 
and 6.8 um and substantial short-wavelength flux with dis- 
tinctive colours. (Such sources could be easily searched for 
in three-colour images consisting of IRAC bands 1, 2, and 
3, where they have a yellow colour in contrast to the blue of 
normal stars and the pink to red of the red sources, where 
bands 1, 2, and 3 are blue, green, and red, respectively.) 
The red sources are dominated by PAHs from 5 to 9 um. (It 
is conceivable that if there should be any 6.0 or 6.8 um ice 
features in the red sources, they would not be detectable be- 
cause of the low signal/noise in the relatively low continuum 
and the contamination by the strong PAH features.) (2) In 
the 9-20 urn range, the outflow sources show substantial 
absorption from both silicates at 10 urn and cold CO2 ice at 
15.2 um, whereas the red sources are extincted only by the 
diffuse ISM, with little CO2 ice absorption. (3) The four red 
sources that are extended and the more luminous all ionize 
their own Hn regions, whereas none of the outflow sources 
produce any ionizing photons (except perhaps in shocks). 
(4) Five of the seven outflow sources produce [S 1] 25.25 u.m 
line emission from their outflow regions. This is indicative of 
low velocity shocks. The two red sources with forbidden-line 
morphology indicating shocks in outflows must have higher 
velocity shocks given the higher excitation of the lines. (5) 
Six of the seven outflow sources produce 6.7 GHz Class II 
methanol maser emission (as does one of the red sources). 
Such masers are commonly found in the earliest stages of 
massive star formation (e.g., Ellingsen 2006; Breen et al. 
2010), and are particularly common among the EGOs stud- 
ied by Cyganowski et al. (2009). 

We suggest that the two groups of YSOs correspond to 
two separate evolutionary stages, with the outflow stage be- 
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ing younger than the red stage. (1) The outflow sources have 
large, thick envelopes containing cold dust grains coated 
with ices. The sizes of the envelopes are some tens of arcsec 
across from the extinction seen in the images (Supporting 
Information Figs. 32, 33, 35 - 38) and the optical depth 
maps, corresponding to a few times 10 4 to > 10 5 au in ra- 
dius. The envelopes of the red sources have smaller optical 
depths since there is little variation in extinction with dis- 
tance from the YSO and the optical depth at the centre is 
not much larger than the optical depths due to intervening 
ISM. (2) The central stars of the red sources have already 
contracted to the point that the protostars' surface temper- 
atures are hot enough to produce photons that can ionize 
hydrogen and neon and doubly ionize sulphur and even some 
Ne ++ . The only ionized gas that could be associated with 
any of the outflow sources is found at the end of the outflow 
of G333.466 and is possibly shock ionized. 

This separation by evolutionary stage is consistent with 
the results of fitting the SEDs with the models of Robitaillc 
et al. (2006) as seen in Table 3 (which indicates that the 
SEDs do have different shapes). For the most part (espe- 
cially for the models with unconstrained interstellar extinc- 
tion), the models of the outflow YSOs have younger ages 
and cooler effective temperatures for the central protostar, 
thereby requiring much larger protostar radii, since the lu- 
minosities and masses have similar ranges and distributions. 
This is in accordance with the models of Hosokawa et al. 
(2010), who find that as a protostar accretes mass, if the 
accretion rate is high enough that the protostar will end 
up as a massive star (~ 10 -3 M yr _1 ), the protostar be- 
comes very large - 'bloated' - and as such is incapable of 
producing ionizing photons. Later, as the star contracts to 
the main sequence, it becomes hot enough to ionize an H n 
region. We conclude that it is likely that the optically thick, 
outflow sources are bloated protostars with large accretion 
envelopes. 

The Australia Telescope Compact Array (ATCA) ob- 
servations of G333. 125—0.562 are consistent with this, find- 
ing weak thermal emission towards the likely protostar 
(G333. 131— 0.560), but no evidence for an Hn region of any 
type (Lo et al. 2011). 

We return to Fig. 1. We did not observe a complete 
sample of all sources in the G333 cloud - there are two 
more 'possible' outflow sources (EGOs) in the region in the 
catalogue of Cyganowski et al. (2008) and a number of red 
sources can be found by looking more closely at the IRAC 
and MIPS images. Cyganowski et al.'s two EGOs and three 
more luminous red sources that are extended at 8 um are 
also plotted in Fig. 1. The plotted sources all lie outside the 
main clusters and Hn regions, which shine in the UV-excited 
PAH emission that dominates IRAC band 4. It is likely that 
the nine tabulated outflow sources are a complete set at 
their level of luminosity, because IRAC band 2 has good 
spatial resolution and is not limited by confusion. One might 
infer from the low inclination angles of the outflow sources 
in Table 3 that there may be many more high-inclination, 
massive, cold YSOs invisible to IRAC or MIPS that should 
be included in our list of possible outflow YSOs. However, 
we do not find such in the Herschel PACS 130-210 urn image 
except possibly two sources confused with known H n regions 
at 24 urn. We conclude from the scarcity of luminous, NIR- 
and MIR-obscured sources that additional work is needed 



modelling this stage of massive star formation. In addition, 
a more extensive comparison of the images and SEDs from 
2MASS and IRAC through PACS and SPIRE wavelengths 
to identify YSOs from their colours will be a fine project for 
future work. 

On the other hand, there are probably more red sources 
within the region of high IRAC band 4 PAH emission - in 
fact, some FIR point sources can be detected in this region 
on the Herschel SPIRE 250 um image. However, the high 
band 4 emission regions are confused and often saturated 
at the longer MIR and FIR wavelengths, making it difficult 
to confirm whether clumpy structure consists of YSOs and 
not compact Hn regions since the same regions also emit 
free-free radiation (Green et al. 1999). We conclude that it 
is not easy to identify YSOs that are past the bloated stage 
and are well on the way to contracting to the main sequence 
but have not contracted so far that they can ionize their own 
H n regions. This probably means that the time-scale for the 
final contraction stage is very short. The red sources that we 
did observe have molecular gas but our limited survey of the 
region including red sources does not find as much molecular 
gas as is found in the outflow sources - this may be another 
indicator of their more evolved nature. 

We also see in Fig. 1 that there are tendencies for the 
YSOs to cluster, with at least three clusters if we include the 
YSO group in G333.125. There can be both outflow and red 
sources in the same cluster, from which we infer that star for- 
mation occurs over a length of time longer than the contrac- 
tion time needed for the smaller red sources like G332.967. 
The outflow sources within the G333 cloud are all associated 
with dense molecular gas (see our previous surveys by Bains 
et al. 2006; Wong et al. 2008; Lo et al. 2009). In these sur- 
veys we have used the Mopra radio telescope to map numer- 
ous molecular transitions including 13 CO, C 18 0, CS, HCN, 
HNC, HCO+, N 2 H+, C 2 H, HC 3 N, CH 3 OH, SO, and SiO, all 
observed with an angular resolution of ~ 35 arcsec. By com- 
paring the tracers of higher density with the lower density 
13 CO tracer, we are able to investigate the spatial and kine- 
matic power spectra of the entire region. While the analysis 
is ongoing, we find from a spatial power-spectrum analysis 
that all tracers have a similar power law distribution, sug- 
gesting well mixed gas with significant turbulent motions 
(M. Cunningham, in preparation). We speculate that the 
observed clusters of YSOs are examples of star formation 
induced by turbulent compression in the outer regions of 
GMCs (e.g., Elmegreen 2007). 



5 SUMMARY AND CONCLUSIONS 

We have mapped 13 YSOs that were identified in the G333 
CMC on Spitzer IRAC and MIPS images with the Spitzer 
IRS from 5-36 pi. We use these spectra plus available pho- 
tometry and images to characterize the YSOs. The objects 
are divided into two groups: YSOs associated with extended 
emission in IRAC band 2 at 4.5 um ('outflow sources') and 
YSOs that have extended emission in all IRAC bands peak- 
ing at the longest wavelengths ('red sources'). 

We find the following results: The source luminosities 
range from a few times 10 3 L to a few times 10 4 L Q . Mod- 
elling of the SEDs indicates that these YSOs have masses 
from 8-25 Mq. The spectral energy distributions (SEDs) 
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of all sources peak between 50 and 110 [im, thereby show- 
ing their young age, but the SEDs of the outflow sources 
peak at slightly longer wavelengths than the SEDs of the 
red sources. 

All spectral maps show ionized forbidden lines and PAH 
emission features (most of this is due to background emission 
from the G333 cloud and must be carefully subtracted from 
the source spectra). For four of the red sources, the line 
emission is concentrated at the centres of the maps, from 
which we infer that these YSOs are the source of ionizing 
photons. The YSO luminosities and low excitation of the 
H II regions are both indicative of early B stars. 

All the YSOs associated with outflows show evidence of 
massive envelopes surrounding the star. The extinction from 
these envelopes is characterized by deep silicate absorption 
features and by absorption by ices at 6.0, 6.8, and 15.2 [im. 
Six out of seven of these YSOs also have 6.7 GHz Class II 
methanol maser emission. At least four of the objects are 
found in clusters of YSOs. Ionized lines are associated with 
two of the clusters containing outflow YSOs; however, the 
gas is more likely ionized by another star in the cluster since 
the lines peak at least 15 arcsec from the outflow YSO. 

For several objects of both types, the lines from the 
most highly excited ions (Ne ++ ) peak at some distance from 
the peak of the low excitation ion lines; we suggest that these 
ions indicate the presence of shocked gas. There is shocked 
gas associated with the 4.5 [im emission, but this is seen 
in the presence of the [S i] line in five of the seven outflow 
sources and not in the presence of any H2 lines, as had been 
suggested by other observers. 

We conclude that this spectroscopic survey has provided 
unique information on this stage of the formation of mas- 
sive stars including the identification of the stars producing 
the IRAC 4.5 u.m outflows as having envelopes cool enough 
to contain ice-coated grains and the characterization of the 
only slightly more evolved YSOs as no longer exhibiting ice 
features but already producing photons energetic enough to 
ionize a small H 11 region. 



ACKNOWLEDGMENTS 

This work is based on observations made with the Spitzer 
Space Telescope, which is operated by the Jet Propulsion 
Laboratory (JPL), California Institute of Technology, under 
a contract with NASA. Support for this work was provided 
by NASA through RSA 1376528 issued by JPL/Caltech. 
The IRS was a collaborative venture between Cornell Uni- 
versity and Ball Aerospace Corporation funded by NASA 
through the JPL and Ames Research Center. SMART was de- 
veloped by the IRS Team at Cornell University and is avail- 
able through the Spitzer Science Center at Caltech. This 
research made use of Tiny Tim/Spitzer, developed by John 
Krist for the Spitzer Science Center. The Center is managed 
by the California Institute of Technology under a contract 
with NASA. This research used observations with AKARI, a 
JAXA project with the participation of ESA. Herschel is an 
ESA space observatory with science instruments provided 
by European-led Principal Investigator consortia and with 
important participation from NASA. We thank S. Molinari 
and the Hi-GAL team for waiving the proprietary period 
such that we could look at their beautiful data in advance 



of publication. We thank J. D. T. Smith for making the 
pahfit program available such that users can modify it for 
other needs. We thank the Spitzer operations team for man- 
aging the telescope schedule such that the cryogens lasted 
until mid-May, 2009. Without this laudable effort, these ob- 
servations would not have been possible. We thank Sean 
Carey and Roberta Palladini for giving us the MIPSGAL 
70 (im images prior to publication and we thank the referee, 
Barbara Whitney, for her thoughtful and thought-provoking 
comments on the manuscript. 



APPENDIX A: LINE FLUX TABLES 

The source fluxes, measured after the subtraction of the 
PAH template, are given in Table Al. The contour levels 
for Figs. 22 - 26 and 28 - 30 are given in Table A2. The 
measured line fluxes needed for Table 5 are given in Table 
A3. 
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Table Al. Measured source continuum fluxes. The continuum fluxes for each source after PAH subtraction 
are given in Jy. G333.466-YSO is the point source at the head of the outflow that has the ice absorption 
features. G333.466-0.164 is the whole source, IRAS 16175-5002. 



Wavelength (urn) 5.48 7.99 9.52 12.02 14.02 18.02 20.01 25.03 30.00 35.02 
Source 



G332 


.725- 


-0.620 a 


0.966 


0.487 


0.810 


2.76 


2.82 


3.88 


6.09 


16.5 


30.7 


55.2 


G332 


.800- 


-0.595 


3.93 


9.76 


14.4 


59.5 


97.5 


200 


293 


526 


709 


946 


G332 


.813- 


-0.700 


0.770 


1.04 


0.635 


2.07 


2.68 


3.36 


4.92 


12.0 


19.8 


31.2 


G332 


.942- 


-0.686 


1.24 


1.66 


0.384 


1.68 


1.50 


1.24 


1.79 


6.91 


17.9 


28.8 


G332 


.963- 


-0.679 


0.897 


1.77 


0.056 


1.09 


2.36 


2.81 


5.98 


36.8 


83.2 


96.5 


G332 


.967- 


-0.683 


0.279 


0.640 


0.241 


1.17 


1.24 


2.18 


3.31 


8.61 


16.3 


24.3 


G333 


.131- 


-0.560 


0.211 


0.264 


0.183 


0.757 


0.832 


1.22 


1.86 


4.59 


8.68 


15.9 


G333 


.163 


-0.100 


0.683 


2.16 


0.235 


5.42 


9.59 


11.8 


20.9 


71.7 


136 


206 


G333 


.184- 


-0.091 


1.13 


1.42 


0.046 


1.46 


2.33 


1.77 


3.32 


12.5 


26.1 


46.3 


G333 


.212- 


-0.105 


0.164 


0.857 


0.088 


1.23 


1.92 


1.92 


3.51 


10.7 


18.0 


27.6 


G333 


.340- 


-0.127 


0.548 


2.81 


0.815 


5.82 


9.85 


17.2 


25.5 


57.1 


82.6 


105 


G333 


.375- 


-0.202 


1.65 


1.69 


0.318 


2.22 


3.33 


5.52 


9.30 


24.9 


40.8 


55.1 


G333 


.466- 


-0.164 


1.84 


4.52 


0.201 


6.26 


11.2 


13.2 


22.2 


78.9 


157 


225 


G333 


.466- 


-YSO 


0.884 


1.74 


0.0039 


1.29 


2.40 


3.02 


5.77 


26.8 


60.1 


88.1 



For G332.725-0.620, the fluxes are measured at 7.40 and 10.18 urn instead of 7.99 and 9.52 urn. 



Table A3. Measured line intensities of the sources producing ionizing photons. The intensities are integrated over a field of view 
comprising all or most of a source (which can be the area emitting forbidden lines in an otherwise outflow source), with the field 
of view, or aperture, tabulated in square arcsec. In addition to the line intensities and errors, the table gives the average optical 
depth at 9.6 urn, to,. 6. 
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Intensities must be multiplied by a factor of 6.3 to compensate for the incomplete sampling. 
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Table A2. Contour levels for Figs. 22 - 26 and 28 - 30. The 
figures noted are contour plots of various line intensities. For each 
figure and line, the contours are usually 0.95, 0.9, 0.8, 0.7, 0.6, 0.5, 
0.4, 0.3, 0.2, 0.1, 0.05, and 0.025 times the maximum observed 
intensity (exceptions are footnoted), where the line, maximum 
intensity, and minimum contour level are given in the table. 
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SUPPORTING INFORMATION 

Additional Supporting Information may be found in the on- 
line version of this article: 

Figures 31-41. In this Supporting Information, we 
show images from the IRAC (panels a - d are IRAC bands 
1-4, respectively) and MIPS (panels e - f are MIPS im- 
ages at 24 and 70 urn images) of each source. An example 
of the Supporting Information figures is shown in Fig. 2 for 
G333. 131—0.560. In each figure, on panels (c) - (e), we over- 
plot the locations of the IRS slits. On panel (b) the locations 
of any 6.7 GHz methanol maser (Caswell 2009) are marked 
with white crosses. All images are shown in logarithmic scal- 
ing. For Fig. 31 of G332.800-0.595, the MIPS 24 urn image 
is replaced by the MSX 21 um image because the MIPS 
24 um image is badly saturated. In Fig. 34, G333. 212-0.105 
is not visible in the MIPS 70 um image (panel f), probably 
because of high background and low Spitzer resolution, but 
it is visible in all Herschel PACS and SPIRE 70 - 600 um 
images. 

This paper has been typeset from a TgX/ WFpJi file prepared 
by the author. 
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Figure 31. IRAC, MSX, and MIPS images of G332.800-0.595. Panel (a): IRAC band 1 (3.6 |im) image. Panel (b): IRAC band 2 
(4.5 |xm) image. Panel (c): IRAC band 3 (5.8 |xm) image with overlays of the SL slit positions. Panel (d): IRAC band 4 (8.0 |xm) image 
with overlays of the SH slit positions. Panel (e): MSX 21 |xm image with overlays of the LH slit positions. The MIPS 24 nm image was 
not used in this figure because it is badly saturated. Panel (f): MIPS 70 um image. 
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Figure 32. IRAC and MIPS images of G332.963-0.679 and G332.967-0.683, which were observed in separate AORs. Panel (a): IRAC 
band 1 (3.6 |xm) image. Panel (b): IRAC band 2 (4.5 urn) image. The white cross marks the location of a 6.7 GHz methanol mascr 
(Caswell 2009). Panel (c): IRAC band 3 (5.8 urn) image with overlays of the SL slit positions. Panel (d): IRAC band 4 (8.0 urn) image 
with overlays of the SH slit positions. Panel (e): MIPS 24 um image with overlays of the LH slit positions. Panel (f): MIPS 70 urn image. 
The centre of the MIPS 24 urn image (panel c) is saturated at the location of G332. 963— 0.679. 
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Figure 33. IRAC and MIPS images of G333. 184—0.091. The images and overlays have the same description as Fig. 32 except that the 
MIPS 24 |xm image is not saturated. Panel (c) also shows the location of the SL1 spectra from AOR 25915648. 
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Figure 34. IRAC and MIPS images of G333. 212— 0.105. The images and overlays have the same description as Fig. 32 except that there 
is no detected methanol maser and the MIPS 24 win image is not saturated. G333. 212— 0.105 is not visible in the MIPS 70 \aa image 
(panel f), probably because of high background and low Spitzer resolution, but it is visible in all Herschel PACS and SPIRE 70 - 600 um 
images. 
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Figure 35. IRAC and MIPS images of G332. 725— 0.620. The images and overlays have the same description as Fig. 32 except that the 
MIPS 24 um image is not saturated. Panel (a): Note that the outflow region is visible on this image as well as on the IRAC band 2 
(4.5 urn) image. 
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Figure 36. IRAC and MIPS images of G332. 813— 0.700. The images and overlays have the same description as Fig 
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Figure 37. IRAC and MIPS images of G332. 942— 0.686. The images and overlays have the same description as Fig. 35. 




Figure 38. IRAC and MIPS images of G333.466-0.164 (IRAS 16175-5002). The images and overlays have the same description as Fig. 
32. The MIPS 24 um image (panel e) is saturated at the location of the maximum of the ionized forbidden line emission, ~ 15 arcsec 
north-west of G333.466-0.164, the YSO at the head of the 4.5 urn outflow (panel b). 
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Figure 39. IRAC and MIPS images of G333. 163— 0.100. The images and overlays have the same description as Fig. 32. Panel (c) also 
shows the location of the SL2 spectra from AOR 25913344. The centre of the MIPS 24 urn image (panel e) is saturated at the location 
of G333.163-0.100. 
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Figure 40. IRAC and MIPS images of G333. 340— 0.127. The images and overlays have the same description as Fig. 32 except that there 
is no detected methanol maser. The centre of the MIPS 24 nm image (panel e) is saturated at the location of G333. 340— 0.127. 
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Figure 41. IRAC and MIPS images of G333. 375— 0.202. The images and overlays have the same description as Fig. 32 except that there 
is no detected methanol maser. The centre of the MIPS 24 um image (panel e) is saturated at the location of G333. 375— 0.202. 



